Olivier Regniers

Lionel Bombrun

Virginie Lafon
email: virginie.lafon@u-bordeaux.fr

Christian Germain

L Regniers

C Bombrun

Supervised Classification of Very High Resolution Optical Images Using Wavelet-Based Textural Features

Keywords: Texture, wavelet, multivariate models, very high resolution, optical image, classification

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

W ITH the launch of more than ten Very High Resolution (VHR) optical satellites in the past 15 years (Quick-Bird, GeoEYE, WorldView, Pléiades, ...), satellite image data of metric and submetric resolution became increasingly available. The level of details provided in such data enables to distinguish geometric structures only observable through their spectral properties at coarser resolutions. These geometric structures can be associated with a regular spatial organisation specific to particular types of land covers. In agricultural landscapes, this is the case of cereal crops, orchards and vineyards which typically display a periodic row structure visible in VHR image data. Because of forest management practices, young tree stands in cultivated forests may also feature specific spatial patterns. Similarly in urban areas, the juxtaposition of buildings can generate a specific framework. These specific spatial patterns can hence be exploited to detect such land covers and improve the classification of VHR optical image data.

Many studies adressed this challenge by considering patterns observed in the landscape as textures. The main objective of texture-based analysis is to explore the local spatial dependencies observed between neighboring pixels in the image. This analysis generally leads to the extraction of a smallsized set of features that can be further used in a classifier. Various approaches were proposed in the literature to represent textures for the classifcation of VHR image data. Among these approaches, the Gray Level Co-occurrence Matrix (GLCM) initially proposed in [START_REF] Haralick | Textural features for image classification[END_REF] is still very popular within the remote sensing community. In many publications, texture descriptors derived from GLCMs were successfully used for various remote sensing applications, e.g. the classification of urban areas [START_REF] Pacifici | A neural network approach using multi-scale textural metrics from very high-resolution panchromatic imagery for urban land-use classification[END_REF], [START_REF] Pesaresi | A robust built-up area presence index by anisotropic rotation-invariant textural measure[END_REF], the mapping of forest species [START_REF] Coburn | A multiscale texture analysis procedure for improved forest stand classification[END_REF], [START_REF] Franklin | Incorporating texture into classification of forest species composition from airborne multispectral images[END_REF], the estimation of forest structure variables in mono-specific forests [START_REF] Béguet | Modelling-based feature selection for classification of forest structure using very high resolution multispectral imagery[END_REF], [START_REF] Champion | Retrieval of forest stand age from SAR image texture for varying distance and orientation values of the gray level co-occurrence matrix[END_REF], [START_REF] Kayitakire | Retrieving forest structure variables based on image texture analysis and IKONOS-2 imagery[END_REF] and the classification of agricultural land covers [START_REF] Warner | Spatial classification of orchards and vineyards with high spatial resolution panchromatic imagery[END_REF], [START_REF] Balaguer | Definition of a comprehensive set of texture semivariogram features and their evaluation for object-oriented image classification[END_REF]. Rather than directly characterizing the texture in the image domain as it is the case with GLCM, other authors suggested to proceed with the texture analysis in a transformed domain of the original data by applying filter banks. For example, texture features extracted by applying scale and orientation selective Gabor filters were proposed in [START_REF] Aksoy | Automatic mapping of linear woody vegetation features in agricultural landscapes using very high resolution imagery[END_REF] to map hedgerows in rural landscapes. An unsupervised segmentation algorithm based on Gabor filters was also introduced in [START_REF] Rabatel | A non-supervised approach using Gabor filters for vine-plot detection in aerial images[END_REF] for the detection of vineyards. In the same way as Gabor filters, wavelet filters also offer a multi-scale and multi-orientation framework for the texture analysis. Features such as energy and entropy [START_REF] Lucieer | Panchromatic wavelet texture features fused with multispectral bands for improved classification of highresolution satellite imagery[END_REF] or GLCM descriptors [START_REF] Ruiz | Texture feature extraction for classification of remote sensing data using wavelet decomposition: a comparative study[END_REF] can be extracted from each wavelet subband to characterize the texture in this transformed domain. In another common approach, probabilistic models are used in the image domain to describe local spatial dependencies and further characterize the textural information. Markov Random Fields, known for their use in the regularization of labeled image, can be modeled with these probabilistic distributions for the classification of VHR remote sensing data [START_REF] Jackson | Adaptive bayesian contextual classification based on Markov random fields[END_REF], [START_REF] Zhao | Classification of high spatial resolution imagery using improved Gaussian Markov randomfield-based texture features[END_REF]. Finally, rather than relying on pre-defined texture features, the increasingly popular deep-learning algorithms efficiently detect patterns in images through unsupervised or semi-supervised feature learning in a deep neural network architecture with many applications in remote sensing data [START_REF] Vaduva | Deep learning in very high resolution remote sensing image information mining communication concept[END_REF], [START_REF] Romero | Unsupervised deep feature extraction for remote sensing image classification[END_REF].

In the past few years, texture features derived from the wavelet domain gained in popularity. The wavelet domain indeed simplifies the way textures can be apprehended as it offers a multi-scale and multi-orientation framework in which each wavelet subband contains a different part of the frequency content of the original image. Low frequencies are discarded and only high frequencies are investigated during the texture analysis. Recently, new approaches emerged to analyze textures by relying on the combined use of probabilistic models and wavelet transform. Probabilistic models used in the wavelet domain can be more simple than in the image domain as unimodal distributions are often associated with wavelet coefficients. This technique consists in modeling the distribution of wavelet coefficients with a probability density function (PDF) whose parameters are estimated according to the maximum likelihood principle [START_REF] Do | Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF]. These parameters are determined for each high frequency subband of the wavelet decomposition and can be further used as a textural signature. Univariate models such as the Generalized Gaussian Distribution [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF] or the Generalized Gamma Distribution [START_REF] Choy | Statistical wavelet subband characterization based on generalized gamma density and its application in texture retrieval[END_REF] were for example proposed to represent the marginal distribution of wavelet coefficients. However, univariate models do not fully describe the statistical behavior of the wavelet coefficients as they do not include informations relative to the spatial dependencies between these coefficients. To achieve this, multivariate models were later introduced to represent the joint distribution of the wavelet coefficients. In this case, all the wavelet coefficients contained in a defined neighborhood within the subband are considered as a local observation of the spatial dependency. All these observations are clustered in a random vector whose distribution can be further modeled using multivariate PDF. The multivariate Gaussian model, the Multivariate Generalized Gaussian Distribution (MGGD) [START_REF] Verdoolaege | Multiscale colour texture retrieval using the geodesic distance between multivariate generalized Gaussian models[END_REF], models based on Spherically Invariant Random Vectors (SIRV) [START_REF] Bombrun | A complete coverage of log-cumulant space in terms of distributions for polarimetric SAR data[END_REF] or copula-based models [START_REF] Kwitt | A joint model of complex wavelet coefficients for texture retrieval[END_REF] are among the most frequently studied models in this context. The classification performances of these multivariate models were mostly investigated in a content-based image retrieval framework applied on natural texture databases such as VisTex [START_REF] Pickard | Vistex database[END_REF] or Brodatz [START_REF] Brodatz | Textures: a photographic album for artists and designers[END_REF] and showed greater retrieval performances than univariate models. Despite the growing interest of the image processing and computer vision communities for this approach, the use of wavelet-based multivariate modeling has not yet been extensively explored for the processing of VHR optical remote sensing data. In a recent work, we demonstrated that waveletbased multivariate models, in particular SIRV-based models, showed promising performances for the classification of forest stand age in mono-specific cultivated forests [START_REF] Regniers | Wavelet-based texture features for the classification of age classes in a maritime pine forest[END_REF].

The potential of the use of multivariate models in the wavelet domain is already known in the image processing and computer vision communities [START_REF] Camps | Chapter 2 -the statistics of remote sensing images[END_REF]. The main contribution of this paper is to demonstrate that such texture analysis approaches are also suitable for the supervised classification of textured soil occupations in VHR optical remote sensing data. Moreover, we propose a complete strategy to apply such models in the context of the classification of VHR optical satellite data. This strategy consists in two steps. First, a content based image retrieval system is used to identify the best probabilistic models to be considered in the context of application. Once identified, the best models are used in a region-wise supervised classification procedure applied on a pre-partitioned image. From a more practical point of view, the main objective of this paper is also to highlight the universality of the proposed strategy which can easily be adapted to various thematic applications with a limited parameters to be set. To this aim, texture-based classification results are presented and discussed for two application examples: the detection of vineyards and the detection of cultivated oyster fields.

The paper is organized as follows. In section II, a more thorough insight on the wavelet transform and the probabilistic modeling of wavelet subbands is given. The proposed strategy for the application of wavelet-based texture features is detailed in section III. This third section is divided in two parts. In the first part, the content-based image retrieval framework is presented with a focus on similarity measurements between multivariate models. The proposed supervised classification procedure is described in the second part. In the following section, texture-based classification results are discussed and compared with a reference approach using GLCM descriptors. Conclusions and perspectives are finally presented in the last section.

II. WAVELET-BASED MULTIVARIATE MODELING

Wavelet transforms provide a multi-resolution analysis framework in which the image is decomposed in a series of subband each containing an orientation and scale specific part of the frequency content of the image. This multi-resolution scheme of analysis is convenient for texture analysis as the pattern constituting the texture can be composed of elements of variable sizes. The most standard wavelet decomposition used for image filtering is the Discrete Wavelet Transform (DWT) initially proposed by Mallat [20]. This decomposition consists in applying combinations of oriented one-dimensional high-pass (H) and low-pass (L) filters to the image resulting in four subbands (three orientation specific high-pass filtered subbands (HH, HL, LH) and one low-pass filtered subband (LL)) subsampled by a factor 2 (dyadic wavelet) at each scale of the decomposition. A subsequent scale is obtained applying the same combinations of filters to the LL subband. The DWT has the advantage of producing independent orthogonal subbands without redundancy of information. Nevertheless, the DWT is limited in two points: the directional selectivity (only three main orientations can be considered) and the shift variance [START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF]. Other configurations of the wavelet decomposition which overcome these issues were later proposed. In wavelet packets [START_REF] Laine | Texture classification by wavelet packet signatures[END_REF], the decomposition is applied at each scale on all the subbands of the previous scales including the high-pass filtered subbands, allowing a finer analysis of the high frequencies in the image. The stationary wavelet transform [START_REF] Unser | Texture classification and segmentation using wavelet frames[END_REF] produces undecimated subbands preventing issues related to shift invariance. Steerable pyramids [START_REF] Simoncelli | The steerable pyramid: A flexible architecture for multi-scale derivative computation[END_REF] enables a shift invariant and multi-orientation decomposition with more than three main orientations. Whichever the chosen wavelet decomposition, the joint distribution of wavelet coefficients in each subband can be further modeled as illustrated in Fig. 1. The aim of this paper is not to compare these various configurations of wavelet decomposition and we mainly focus on the comparison of performances of different multivariate models applied on DWT subbands.

To conduct the modeling of the spatial dependencies within a subband, the wavelet coefficients contained in the neighborhood of the current spatial position are clustered in a random vector k. Several multivariate models were proposed to represent the distribution of the random vector k. We first focus on the well-known multivariate Gaussian model and further introduce two generalizations of this model: the SIRV-based models and the copula-based models. For each of these models, the technique used to estimate the distribution parameters is detailed in the following paragraphs. These estimated parameters η are further used as a textural signature in the classification procedure.

A. Multivariate Gaussian Model

If we assume the gaussianity of the dependency to be modeled, the observed vector k is a real Gaussian vector of zero mean and covariance matrix [M]. A zero mean can be considered here as only the high frequency subbands are modeled. This stochastic Gaussian distribution is then characterized by an estimation of the covariance matrix [M] using the Sample Covariance Matrix estimator (SCM):

[ M] SCM = 1 N N i=1 k i k i ( 1
)
where N is the number of observations and k i is the transposed vector of k i . Hence, the textural signature obtained from this model contains an estimation of the covariance matrix for each subband of the decomposition. However, it is well known that Gaussian models are not flexible enough to capture the "heavy-tail" behaviour generally observed in the distribution of wavelet coefficients [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]. Other multivariate models offering more adjustement flexibility were therefore proposed as generalizations of the Gaussian case. The SIRV-based and copulabased models are among such generalizations.

B. SIRV-based Models

In SIRV-based multivariate models, the generalization of the Gaussian case is obtained as follows. The observed vector k is decomposed as k = √ τ z where z is a Gaussian vector of zero mean and covariance matrix [M] and τ a positive random scalar variable called the multiplier, τ and z being assumed independent. The multiplier τ enables to modulate the distribution of the Gaussian vector z to better fit the heavytail behaviour of the distribution of wavelet coefficients. The PDF of the random variable τ is not specified in the definition of SIRV models and various univariate models derived from the Pearson system of distributions were proposed to represent the multiplier [START_REF] Bombrun | A complete coverage of log-cumulant space in terms of distributions for polarimetric SAR data[END_REF]. In this paper, we will focus mainly on two particular configurations of the SIRV models. First, when considering a deterministic multiplier τ , a Gaussian-like distribution is produced (SIRVgauss). Here, the multiplier τ varies from one sample k i to another and can be represented by a Dirac [START_REF] Pascal | Performance analysis of covariance matrix estimates in impulsive noise[END_REF]. Second, for an Inverse Gamma distributed multiplier τ , the SIRV vector k follows a G 0 distribution (SIRVg0) [START_REF] Freitas | The polarimetric G distribution for SAR data analysis[END_REF].

The SIRV representation is not unique and a normalization condition is necessary. Classicaly, the normalization is imposed on the trace of the covariance matrix, set here to p the dimension of the target vector k. An approximated maximum likelihood estimator of the normalized covariance matrix is the solution of the following recursive equation known as the "Fixed Point" estimator (FP) [START_REF] Vasile | Coherency matrix estimation of heterogeneous clutter in high-resolution polarimetric sar images[END_REF]:

[ M] F P = f ([ M] F P ) = p N N i=1 k i k i k i [ M] -1 F P k i with T r([ M] F P ) = p (2)
The existence and the uniqueness of the Fixed Point estimator as well as the convergence of the recursive algorithm whatever the initialization was established in [START_REF] Pascal | Covariance structure maximum-likelihood estimates in compound gaussian noise: Existence and algorithm analysis[END_REF]. When used in SIRV configurations, the FP estimator has been proven to be more robust than the SCM estimator used for the multivariate Gaussian model [START_REF] Pascal | Performance analysis of covariance matrix estimates in impulsive noise[END_REF]. For a given covariance matrix [M], the maximum likelihood estimator of the multiplier τ i is then defined by:

τi = k i [M] -1 k i p (3)
In both SIRVgauss and SIRVg0 cases, the obtained textural signature contains the FP estimator of the covariance matrix, completed in the case of the SIRVg0 model with the parameters of the univariate Inverse Gamma PDF chosen to represent τ (one shape and one scale parameter).

C. Copula-based Models

Another generalization of the Gaussian case is given by the family of copula-based models. Here, the margins of the observed p-dimensional vectors k = [k 1 , . . . , k p ] are first modeled independently by a chosen univariate PDF. A transformation is then applied on the modeled marginals by computing a cumulative distribution function CDF F 1 , . . . , F p resulting in a set of transformed variables (F 1 (k 1 ), . . . , F p (k p )). Then, a copula C is used to represent the relationship between these transformed variables. The definition of the copula does not specify the distribution associated with the CDF or the marginals and various configurations of copulas are possible [START_REF] Kwitt | A joint model of complex wavelet coefficients for texture retrieval[END_REF], [START_REF] Stitou | Copulas based multivariate gamma modeling for texture classification[END_REF], [START_REF] Lasmar | Gaussian copula multivariate modeling for texture image retrieval using wavelet transforms[END_REF]. The advantage of copulas in comparison to SIRV models is that they enable to define a joint distribution with arbitrary marginals.

The definition of the copulas is derived from the Sklar's theorem [START_REF] Sklar | Random variables, joint distribution functions, and copulas[END_REF] which states that given a p-dimensional random vector k = [k 1 , . . . , k p ] with a CDF F and continuous marginals F 1 , . . . , F p , there exists a unique copula C such that:

F (k 1 , . . . , k p ) = C(F 1 (k 1 ), . . . , F p (k p )) (4)
Given that the CDF F and the copula C are continuous, the joint density f (k) can be expressed as follows:

f (k 1 , . . . , k p ) = c(F 1 (k 1 ), . . . , F p (k p )) p r=1 f r (k r ) ( 5
)
where c is the copula density and f r , r = 1, . . . , p are the marginal PDFs. Among the configurations of copulas, we will focus on the use of a multivariate Gamma distribution with a Gaussian copula (GGC) as proposed in [START_REF] Stitou | Copulas based multivariate gamma modeling for texture classification[END_REF]. In this configuration, the copula density c is modeled by a Gaussian distribution and the marginals f are modeled with univariate Gamma distribution. Similarly to SIRV models, the multivariate Gaussian part of the copula model is modulated by the univariate distribution applied to the marginals in order to better capture the non-Gaussianity of the distribution of the wavelet coefficients.

In the case of Gaussian copulas, the maximum likelihood estimation of the parameters can be done in two steps. The first consists in the estimation of the parameters η = (η 1 , . . . , ηp ) of each marginal f 1 , . . . , f p independently:

ηr = arg max ηr log p r=1 f r (k r,i , η r ) (6)
The estimation of the covariance matrix [M] is then obtained in a second step after transforming the observed vectors k i into Gaussian vectors k G i . The maximum likelihood estimator of the matrix [M] is then defined as:

[ M] = 1 N N i=1 k G i k G i (7)
Hence, the extracted textural signature obtained with Gaussian copula models is also composed of an estimation of a covariance matrix and the parameters of the univariate distribution chosen to represent the marginals (one shape and one scale parameter in the case of the Gamma distribution). Further details on the estimation of these parameters can be found in [START_REF] Stitou | Copulas based multivariate gamma modeling for texture classification[END_REF], [START_REF] Lasmar | Gaussian copula multivariate modeling for texture image retrieval using wavelet transforms[END_REF]. The four multivariate models described above (SCM, SIRVg0, SIRVgauss, GGC) were tested in two different contexts of VHR data classification. Before presenting these results in section IV, the strategy proposed to include these wavelet-based multivariate models in a supervised classification procedure is explained in the next section.

III. METHODOLOGY

Prior to the supervised classification, a first analysis using a Content-Based Image Retrieval framework [START_REF] Smeulders | Content-based image retrieval at the end of the early years[END_REF] is conducted Fig. 2: Architecture of the Content-Based Image Retrieval framework (adapted from [START_REF] Do | Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF]). on a pre-existing texture database in order to identify which are the multivariate models capable to most efficiently distinguish the considered texture classes. The concept of distance or similarity measurement between multivariate models is inherent to the CBIR framework and is defined here below. Once identified, the best models are further used in a supervised classification procedure described in the second part of this section.

A. Content-Based Image Retrieval

The Content-Based Image Retrieval framework (CBIR) is a universal technique aiming at the retrieval of images from a database depending on color, shape or texture features describing the visual content of these images. In this study, we focus on the use of textural features derived from wavelet-based multivariate models. This general framework is frequently used in the literature to compare the ability of probabilistic models to represent the textural content of image patches belonging to the VisTex or Brodatz databases [START_REF] Pickard | Vistex database[END_REF], [START_REF] Brodatz | Textures: a photographic album for artists and designers[END_REF]. In order to apply such framework to a real practical case of image classification, a database of image patches is first created by extracting texture-wise homogeneous regions representative of the different classes of interest from the VHR optical data. Two major tasks are considered in the CBIR: the feature extraction and the similarity measurement (Fig. 2).

During the feature extraction step, the textural signature of all the patches of the database are computed according to the chosen texture analysis approach. In the case of the GLCM for example, a signature consists of a vector of descriptors derived from the co-occurrence matrix computed on one patch, whereas in the case of multivariate models, the signature contains the estimated paramaters of the modeled distribution obtained from each of the wavelet subband as described in the previous section. Once the signatures are computed, one patch is extracted from the database and considered as a query image. A similarity measurement is then calculated between the signature of the query image and all the signatures of the rest of the database, allowing to retrieve the top closest images. A textural analysis method is considered more efficient than another if amongst the top closest images, there is a higher proportion of patches belonging to the same class as the query image. This operation is repeated by considering successively each patch of the database as a query image. In the case of models derived from wavelet subbands, it is worth noting that, by assuming the hypothesis of decorrelation between wavelet subbands, the similarity measured between two images I 1 and I 2 is equal to the sum of similarities measured between the models derived from pairs of subbands W 1 and W 2 at the same scale sc and the same orientation or of the decomposition:

D(I 1 I 2 ) = N sc sc=1 N or or=1 D(W 1sc,or W 2sc,or ) (8)
where N sc and N or are respectively the number of scales and the number of orientations of the decomposition. When using for example textural features derived from a GLCM, the notion of similarity between two images comes down to the measure of a distance between two vectors in the parameter space. A Euclidean or a Mahalanobis distance can be used in such case. For applications based on probabilistic models, Euclidean or Mahalanobis distances are not appropriate as they do not take advantage of the probabilistic nature of the parameters. The use of the Kullback-Leibler Divergence (KLD) has been widely advocated in this context [START_REF] Do | Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF]. The KLD between two models P (k; η 1 ) and P (k; η 2 ) is indeed defined according to the chosen PDF as follows:

KLD = D(P (k; η 1 ) P k; η 2 )) = P (k; η 1 ) log P (k; η 1 ) P (k; η 2 ) dk (9)
where η refers to the model parameters. Unfortunately, a closed form of the KLD is not currently available for all multivariate models. The use of the Rao Geodesic Distance (GD) was proposed as an alternative to the KLD (e.g. [START_REF] Verdoolaege | Multiscale colour texture retrieval using the geodesic distance between multivariate generalized Gaussian models[END_REF], [START_REF] Bombrun | Multivariate texture retrieval using the geodesic distance between elliptically distributed random variables[END_REF]). This distance relies on the Fisher information matrix and is defined as the length of the shortest path connecting two points in the space of the parameters of the distribution. The use of the GD is also preferred since unlike the KLD, it satisfies the triangular inequality. Besides, it was shown that the GD could lead to higher classification performances than the KLD [START_REF] Verdoolaege | Multiscale colour texture retrieval using the geodesic distance between multivariate generalized Gaussian models[END_REF], [START_REF] Lenglet | Statistics on the manifold of multivariate normal distributions: Theory and application to diffusion tensor MRI processing[END_REF]. However, in the case of the GGC model, only a closed form of the KLD is currently available and is used to measure similarities for this model. The approach used to estimate similarities are described here below for the four considered multivariate models. 1) Similarity between Gaussian models: For Gaussian distributions, models are only characterized by their covariance matrix estimated in this study either with the SCM estimator for the multivariate Gaussian model (SCM) or by the FP estimator for the Gaussian SIRV model (SIRVgauss). The GD between two Gaussian models defined by the matrices [M] 1 and [M] 2 of dimension p × p is expressed as follows:

GD([M] 1 [M] 2 ) = p i=1 (ln λ i ) 2 (10)
where λ i are the eigenvalues of

[M] -1 1 [M] 2 .
2) Similarity between SIRVg0 models: Similarity measurements between two SIRVg0 models depends not only on the covariance matrix but also on the other parameters of the distribution. Various strategies were suggested to approximate the GD in this context. Here, we applied the strategy proposed by [START_REF] Bombrun | Multivariate texture retrieval using the geodesic distance between elliptically distributed random variables[END_REF] in which elements of the Fisher information matrix are expressed only as a function of the two shape and scale parameters of the Inverse Gamma distribution. The GD is then approximated by assuming the geodesic coordinate functions as straight lines. Another strategy considering fixed shape parameters was also introduced in [START_REF] Bombrun | Multivariate texture retrieval using the geodesic distance between elliptically distributed random variables[END_REF] for SIRVg0 models, but it was not used in this study as it was shown to be less efficient.

3) Similarity between GGC models: For copula-based models, there are currently no available closed form of the GD. A Monte-Carlo approach to approximate the KLD between two copula-based models was suggested in [START_REF] Kwitt | A joint model of complex wavelet coefficients for texture retrieval[END_REF]. But, this technique often implies high computational costs. More recently, a closed form of the KLD was introduced for Gaussian copulas [START_REF] Lasmar | Gaussian copula multivariate modeling for texture image retrieval using wavelet transforms[END_REF]. This expression equals to the sum of one KLD measured between two Gaussian models for the copula and p KLD between univariate models for the marginals. The symetric version of this KLD was used here for comparing two copulabased models.

In order to assess the retrieval performances in CBIR, various indices can be used (see [START_REF] Müller | Performance evaluation in content-based image retrieval: overview and proposals[END_REF] for a review). Between these, the most commonly used indices are Precision and Recall, usually presented as a Precision versus Recall graph (PR graph). The Precision corresponds to the average fraction of retrieved patches that are relevant amongst the n top closest patches, whereas the Recall is the average fraction of relevant patches that are retrieved in the top n closest image. The PR graph is constructed by varying the parameter n between 1 and N T , the total number of patches in the database. This is equivalent to varying the distance threshold between the query image and the other images of the database. The Recall increases with n, whereas the Precision decreases when n increases. A texture model is therefore considered more performant than others when its associated PR curve is closer to the upper-right corner of the graph. This graphic representation is similar to ROC curves which plots the true positive rate versus the true negative rate. Nonetheless, as recommended in [START_REF] Davis | The relationship between Precision-Recall and ROC curves[END_REF], the use of PR graphs was preferred as it can be more informative when the analysed database is unbalanced in terms of number of samples per class. The PR graph gives a general insight on the ability of a model to detect all the images of a same class without having to scan the entire database. The interpretation of PR graphs is used here to identify which multivariate model would be the best choice for the classification of VHR remote sensing data. In the next part of this section, the proposed strategy for using wavelet-based multivariate models in a supervised classification is presented.

B. Supervised Classification

The classification of VHR optical data using texture features can be addressed in different manners. The most intuitive way is to proceed with a pixel-wise approach using a sliding window to characterize local spatial dependencies. This approach is often preferred when using GLCM descriptors resulting in images of parameters subsequently used as inputs in a classifier [START_REF] Béguet | Modelling-based feature selection for classification of forest structure using very high resolution multispectral imagery[END_REF], [START_REF] Warner | Spatial classification of orchards and vineyards with high spatial resolution panchromatic imagery[END_REF]. This classification technique may have however some drawbacks: (1) the optimal size of the sliding window is not easily determined automatically and trials and errors are often needed to fix it; (2) the VHR of optical data readily translates into high local variability which can be difficult to fully capture during the learning phase without requiring a large database; (3) the computational cost of a pixel-wise analysis can be prohibitive given the size of the objects to classify in regard to the spatial resolution of the image; (4) the use of a sliding window can be the source of classification errors especially at the vicinity of object edges where mixed textures may be present. For all these reasons, we opted for a region-wise classification based on a presegmentation of the image to classify as advocated for example in [START_REF] Gueguen | Classifying compound structures in satellite images: A compressed representation for fast queries[END_REF], [START_REF] Fang | Spectral-spatial classification of hyperspectral images with a superpixel-based discriminative sparse model[END_REF]. This alternative enables to reduce the number of entities to classify to the number of regions obtained by the segmentation and thereby to reduce the computational cost and to increase the homogeneity of the classification results. The proposed supervised classification method operates as follows.

During the learning phase, texture features are extracted from each patch of a learning texture database. This database is the same as the one used during the CBIR analysis. As a reminder, this database was created by extracting texture-wise homogeneous regions representative of the different classes of interest. In parallel, the same texture features are extracted from each region identified during the pre-segmentation. In the case of multivariate models, the wavelet decomposition is applied on the entire image to classify. At each scale of the decomposition, a label image derived from the presegmentation is downsampled to the size of the wavelet subbands so that the wavelet coefficients corresponding to each region can be easily identified. The multivariate modeling is next performed by only considering the neighborhoods contained in each region within each subband. Once texture features are extracted from the learning database and form the regions of the pre-segmentation, a classifier is used to determine which class is affected to each region according to its decision rules. The classification results is further validated by comparing it with a ground truth if available.

Any kind of classification algorithms can be adapted to this classification scheme as long as it includes a notion of similarity or if it respects the probabilistic nature of the proposed multivariate models. Three types of classification algorithms are tested in the following section: k-Nearest Neighbors (k-NN), maximum likelihood (ML) and Support Vector Machine (SVM). The k-NN relies only on a similarity measure between the texture features derived from the learning samples and the region to classify. The majority class among the k-NN is affected to the region. Similarly to the CBIR analysis, a GD or KLD is used to measure similarities. In the ML classifier, a log-likelihood criterion is computed between the multivariate models estimated on the learning samples and the local observations of the spatial dependency within each region. The model learned on the regions is not used in this case. A region is classified according to the class of the sample that maximizes this likelihood criterion. The advantage of the ML classification in comparison to k-NN is that it does not require to re-estimate new features from each of the region. This can be helpful in the case of small-sized regions (over-segmentation) for which the number of local observations in each region may not be sufficient to estimate a stable model. Nonetheless, the computational cost of the ML classifier can be significantly higher, especially when a large number of regions are to be classified. To use a SVM classifier, a nonlinear kernel-based transformation is first applied on the data to project them in a new space where a hyperplan between classes can be defined. To do so, a Gaussian kernel K is adapted to the similarity measure between models by replacing the Euclidean distance with the KLD or the GD [START_REF] Vasconcelos | The Kullback-Leibler kernel as a framework for discriminant and localized representations for visual recognition[END_REF], [START_REF] Jayasumana | Kernel methods on the riemannian manifold of symmetric positive definite matrices[END_REF]:

K(x, x ) = exp - distance(x, x ) 2 2σ 2 (11)
where (x, x ) is a couple of observations and σ the standard deviation of the Gaussian kernel. The value of σ is fixed here to maximise the performance of the classifier through a crossvalidation procedure.

IV. EXAMPLES OF APPLICATIONS

In this section, we illustrate the proposed strategy with two examples of application of wavelet-based multivariate modeling for the classification of VHR optical remote sensing data. Both examples aim at the detection of one particular land cover in VHR panchromatic Pléiades data: vineyards and cultivated oyster fields.

A. Detection of Vineyards

Because of the periodic aspect of vine rows, textural analysis for the detection of vineyards is often performed through the use of frequency-based methods relying on the Fourier transform [START_REF] Wassenaar | Vineyard identification and description of spatial crop structure by per-field frequency analysis[END_REF], [START_REF] Ranchin | An automatic method for vine detection in airborne imagery using the wavelet transform and multiresolution analysis[END_REF], [START_REF] Aksoy | Automatic detection and segmentation of orchards using very high resolution imagery[END_REF]. Among these frequencybased methods, Rabatel et al. [START_REF] Rabatel | A non-supervised approach using Gabor filters for vine-plot detection in aerial images[END_REF] suggested an unsupervised segmentation approach in which each vine plot is isolated by selecting its frequency response in the Fourier domain using adapted Gabor filters. To apply frequency-based approaches in this context, a good contrast is needed between vine rows and inter-rows for the periodic structure to be visible. But, according to the Shannon-Nyquist theorem, such structure would only be observable if the spatial resolution of the image is at least twice as small as the period of the pattern. In many French wine-growing regions and in particular in the area of Bordeaux, the plant density is such that the inter-row distance can be close to 1 meter. Using VHR panchromatic data such as those produced by the Pléiades sensor (spatial resolution of 0.7 m) can therefore be problematic. Aliasing artifacts can be common in these data. This will result in a significantly altered texture in which the periodic aspect of vine rows is no longer clearly visible. In this application, we explore the potentialities of four wavelet-based multivariate models (SCM, SIRVgauss, SIRVg0, GGC) to detect vineyards despite the presence of these artifacts. The classification results obtained with these models are compared with those produced when using GLCM descriptors and the unsupervised segmentation approach proposed in [START_REF] Rabatel | A non-supervised approach using Gabor filters for vine-plot detection in aerial images[END_REF].

To achieve this, we used four panchromatic images acquired by the Pléiades sensor on three wine appellations within the wine-growing region of Bordeaux: Sauternes, Saint-Émilion and Pessac-Léognan. Plant densities in these three appellations are such that the inter-row distances (IRD) are of approximately 1 to 1.2 m in Pessac-Léognan and 1.4 to 2 m in Saint-Émilion and Sauternes. Besides, if the spatial resolution of panchromatic Pléiades data is of 0.7 m, the Ground Sampling Distance (GSD) at the time of the image acquisition only reaches this resolution at the nadir of the satellite for an incidence angle close to 0°. For the Pléiades data used in this study, the incidence angle is significantly higher (from 13°to 33°) except for the Sauternes image (4°) and the GSD is hence significantly increased (from 0.75 to 0.9 m). Aliasing effects and the resulting alteration of textures is all the more important that the GSD is high and that the IRD is low.

From each of the four panchromatic data, a texture database is created by extracting image patches (128×128 pixels) from texture-wise homogeneous areas in the data. These selected patches are representative of the different land cover classes to distinguish in the landscape. In the wine-growing region of Bordeaux, the landscape is dominated by vineyards standing alongside typical peri-urban land covers. Four main classes are therefore considered: vineyards, forest, low vegetation (i.e. meadows, regenerating vineyards, bare soil, etc.) and built areas (Fig. 3). The texture databases constructed from each of the four panchromatic data contains 25 patches for the forest, low vegetation and built areas classes. For the vineyards, the lack of rotation invariance of the wavelet-based texture features has to be taken into consideration and a fairly exhaustive set of patches in terms of orientation has to be created. To do so, we included in the database 25 image patches of vineyards for each main classes of orientation defined by intervals of 20°between 0°and 180°(9 classes of orientation). Each database contains a total of 300 image patches with 225 for vineyards only.

The created databases are first exploited in a CBIR analysis to compare the retrieval performances of the tested multivariate models with those obtained with GLCM descriptors. In order to extract GLCM descriptors from image patches, the distance between pairs of pixels is set to 1 and four Haralick descriptors are selected: entropy, homogeneity, correlation and mean. These four descriptors are extracted from each GLCM computed in four orientations (0°, 45°, 90°, 135°) and further averaged to obtain rotation invariant descriptors [START_REF] Haralick | Statistical and structural approaches to texture[END_REF]. These parameters are selected as they maximize the retrieval performances of the GLCM. Besides, multivariate models are applied on subbands obtained with a DWT with Daubechies db4 wavelet filters. Wavelet db4 filters are chosen as they enable a compromise between the spatial extent of the filter and its localization in frequency. The number of scales in the decomposition is set to 2 producing six wavelet subbands (2 scales × 3 orientations). Neighborhoods of 3×3 are considered within each subband to represent the observations of the local spatial dependency. A total of six sets of model parameters are extracted for each image patch of the database (one set for each wavelet subband).

The PR curves of the wavelet-based multivariate models and the GLCM descriptors are shown in Fig. 4 for each of the texture databases derived from the four Pléiades panchromatic data. It appears clearly on these curves that the retrieval performances of the multivariate models depend directly on the ratio GSD/IRD associated with the processed image data. The more this ratio increases, the more the PR curves struggle to get closer to the upper right corner of the graph, synonym of a perfect retrieval behavior. Indeed, when this ratio is lower than 0.5 (Fig. 4a), i.e. when the conditions of the Shannon-Nyquist theorem are respected, the texture associated with vineyards is less altered and retrieval performances are high with PR curves close to the upper right corner of the graph. In this case, GLCM descriptors as well as models like SCM and GGC are the more performant. Conversely, when the ratio is close to 0.5 (Fig. 4b) or slightly higher than 0.5 (Fig. 4c), the texture is potentially more altered and retrieval performances decrease. In those cases, multivariate models such as SCM and GGC display higher performances than GLCM descriptors. Finally, when the ratio GSD/IRD is significantly higher than 0.5 (Fig. 4d), the level of alteration of textures is such that none of the tested texture analysis succeeds to retrieve the textures efficiently. This lower sensitivity of wavelet-based multivariate models to the texture alteration confirms their ability to capture more informative components of the texture in the wavelet domain than GLCM descriptors in the image domain. This is in agreement with previous results [START_REF] Ouma | Urban features recognition and extraction from very-high resolution multi-spectral satellite imagery: A micro-macro texture determination and integration framework[END_REF] in which waveletbased features were shown to be more adapted to capture lower and mid-frequency texture information, while the GLCM is more efficient to describe high-frequency textural components. Two multivariate models display significantly better retrieval performances than the two others in this particular context of application, i.e. the SCM and the GGC models. These two models are further tested in the proposed supervised classification strategy.

The texture patches database used previously in the CBIR analysis are further used in the classification process as a training database. Prior to the classification, a pre-segmentation is operated on the image using a Mean-Shift algorithm [START_REF] Comaniciu | Mean shift: A robust approach toward feature space analysis[END_REF]. As this algorithm relies directly on the distribution of gray levels, a smoothing filter is first applied on the image to prevent regions from following vine rows, creating thereby irregular narrow contours. The minimum region size parameter of the Mean-Shift algorithm is also set at a low value (500 pixels) so that regions containing mixed textures are limited. The result can be therefore over-segmented. Note that other segmentation algorithms could be used for this step. The identified regions are then classified according to the chosen classifier by either extracting new features from these regions when using k-NN and SVM or by estimating a likelihood criterion when using the maximum likelihood (ML) classifier. Classification results are then compared with a ground truth for validation. The classification results obtained with GLCM descriptors as well as with the SCM and GGC wavelet-based multivariate models are presented in Fig. 5 for an image extract of one of the Pléiades data acquired on the Pessac-Léognan appellation with a GSD/IRD ratio equal to 0.61 (Fig. 5a). These results are compared with the segmentation result achieved with the unsupervised method based on the use of Gabor filters [START_REF] Rabatel | A non-supervised approach using Gabor filters for vine-plot detection in aerial images[END_REF]. Regions classified or segmented as vineyards are displayed in accordance with their consistency with the ground truth (Fig. 5c) by describing them in terms of true positives, false positives, true negatives and false negatives. Only the results obtained with the most efficient classifier are presented for each textural analysis appoach in Fig. 5. The corresponding accuracy indices as well as computational time are shown in Table I. These results are discussed here below. 1 These durations were obtained with an Intel Core i5-3570 CPU at 3.40 GHz and 8 GB of RAM, running the 64 bit version of Windows 7 operating system. The codes for the supervised classification (SVM and ML) and the Gabor filters approach were implemented and run in MATLAB (version 8.5, R2015a, 64 bit). These durations do not include the Mean-Shift presegmentation and the model training, these steps being identical for all texturebased approaches.

Methods using texture features (GLCM and multivariate models) reveal higher detection performances than the segmentation method using Gabor filters. This lesser ability of Gabor filters to detect vine plots in this context is mainly due to the short IRD (1 to 1.2 m) which significantly contributes to the attenuation of the frequency response of the vine plots. Plots with the lowest IRD are practically undetected by this approach (see false negatives in Fig. 5d). Besides, the presence of missing vines in most of the plots interrupts locally the frequency patterns. This translates in small false negative areas appearing in the detected plots. Nevertheless, despite this under-detection problem of Gabor filters, the confusion with other land covers (false positives and user's accuracy) is low in comparison to the use of texture features. Confusions however occur with the neighboring road network which also produces a frequency response in the Fourier spectrum. Concerning the textural approaches, higher detection performances are observed when using multivariate models. Results obtained with these models (Fig. 5f-5g) display indeed higher true positive and lower false positive rates than those obtained with GLCM descriptors (Fig. 5e). This confirms the higher sensitivity to texture attenuation observed for the GLCM in the CBIR analysis. The best compromise between higher true positive and lower false positive is achieved when using the SCM model (Fig. 5f). Most of the false positives are due to confusions with low vegetation which can share similarities with vineyards especially when the texture pattern associated with vine rows is significantly altered. Also, regions defined by the pre-segmentation algorithm tend to spread beyond the limits of vine plots and create mixed regions, especially around narrow paths delimiting the plots. These detection errors are due to the Mean-Shift algorithm rather than to the texture analysis. Overall, wavelet-based multivariate models and the proposed strategy to apply them in this particular context of vineyards classification proved to produce higher detection performances than standard approaches such as GLCM descriptors or Gabor filters. These observations are confirmed by the accuracy evaluation (Table I) with higher accuracy indicators for multivariate models than for GLCM descriptors or Gabor filters. It can however be noted that the computational duration of the classification is, as expected, significantly higher for wavelet-based mulvariate modeling with a computational cost up to 40 times higher than the use of GLCM descriptors and up to 200 times higher than the use of Gabor filters. This higher computational time is mainly due to the calculation of the maximum likelihood criterion between the multivariate models estimated on each training sample and the local observations within each region to classify.

In the following application, the proposed classification strategy is adapted to the detection of small and narrow structures in the case of cultivated oyster fields mapping.

B. Detection of Cultivated Oyster Fields

In the Arcachon Bay, a mesotidal lagoon located in the South-West of France, oysters are traditionally grown on racks of a few tens of meters to more than hundred meters long. These racks are set directly on the sediment in groups of two to five spaced by 1 to 2 meters. These groups of parallel structures are directly visible in VHR optical data and can be detected by using texture analysis in order to assess the evolution of the area occupied by this shell farming industry. Some recent studies demonstrate the interest of the international community to detect and map oyster habitat in tidal flats [START_REF] Choe | Detection of oyster habitat in tidal flats using multi-frequency polarimetric SAR data[END_REF], [START_REF] Gade | Multi-frequency SAR data help improving the monitoring of intertidal flats on the German North Sea coast[END_REF]. In these studies, oyster reefs are mainly detected by making use of the radar backscattering properties specific to these land covers in SAR remote sensing data. Here, we adapt the proposed classification strategy using waveletbased multivariate models to detect such structures in VHR optical data.

For this application, we used a panchromatic image acquired by the VHR Pléiades sensor in April 2013. This image covers a large central part of the Arcachon Bay where the oyster farming activity is predominant. The learning database used for this application is created differently than in the previous application. The use of same-sized image patches to learn texture features may indeed be less appropriate in the case of oyster fields mapping as narrow groups of oyster racks may be separated by a few meters large spaces of bare sediment. The textural signature derived from such patches would therefore include textural informations relative to both oyster racks and bare sediment. Instead, unevenly-sized Regions Of Interest (ROI) whose contours are adapted to the shapes of the structures to detect are extracted from the image. Three classes of land covers are considered: cultivated oyster fields (40 ROIs), abandoned fields (20 ROIs) and foreshores (20 ROIs) (Fig. 6). The latter corresponds to any area of bare sediment where no oyster is grown. For this class, image patches of the same size are extracted since the texture associated with foreshores is homogeneous in vast areas.

First, the learning database is analysed in a CBIR framework to compare retrieval performances and identify which multi- To extract GLCM texture descriptors from ROIs, the same parametrization of the GLCM than in the previous application is used, i.e. the distance between pairs of pixels is set to 1 and four Haralick descriptors (entropy, homogeneity, correlation and mean) are computed in four orientations before being averaged. Besides, multivariate modeling is applied on the six subbands obtained by applying a DWT with Daubechies db4 wavelet filters and 2 scales and 3 orientations of decomposition. Neighborhoods of 3×3 are considered within each subband to represent the observations of the local spatial dependency. As in the previous application, a total of six sets of model parameters are extracted for each image patch of the database (one set for each wavelet subband).

The computed PR curves show higher retrieval performances for almost all the multivariate models in comparison to GLCM descriptors with the exception of the SIRVgauss model (Fig. 7). It is worth noting that conversely to the application dedicated to the vineyards detection, PR curves are decreasing faster and are more distant from the upper right corner. This gives hints on the potentially higher level of confusion between the considered classes. This could be explained by the shared anisotropy between the cultivated oyster fields and abandoned fields which could be the source of classification errors. The performances of the three best multivariate models (SCM, SIRVg0 and GGC) are further assessed in classification.

The ROI database used previously in the CBIR analysis is further used in the classification process as a training set. Prior to the classification, a Mean-Shift pre-segmentation is performed on the image to classify. As groups of oystergrowing racks can be narrow, the Mean-Shift algorithm is deliberately parametrized to create an over-segmentation with the minimum region size parameter set to 100 pixels. Small regions are hence produced (Fig. 8b). The identified regions are then classified according to the chosen classifier by either extracting new features from these regions when using k-NN classifier for GLCM descriptors or by estimating a likelihood The produced classification results show an equivalent capacity of two of the selected multivariate models (SIRVg0, GGC) to correctly detect oyster racks with a low proportion of false positives (Fig. 8f-8g). Conversely, the result obtained with the SCM model (Fig. 8e) reveals higher confusions with the bare sediment spaces between oyster racks with significantly higher false positives in those areas. However, as the two other multivariate models (SIRVg0 and GGC), it prevents confusion errors with abandoned fields. On the other hand, even though it achieves a satisfactory delimitation of cultivated oyster fields without significant confusion with bare sediments, the use of GLCM descriptors also induces classification errors with abandoned fields as it can be seen in the upper right part of Fig. 8d. As it was observed for the detection of vineyards, the best compromise between high level of true positives and low level of false positives is observed in this application when using wavelet-based multivariate models. This also confirms the better performances achieved with these models in the CBIR analysis in comparison to GLCM descriptors. These observations are confirmed by the accuracy evaluation (Table I) with higher accuracy indicators for the SIRVg0 and GGC models than for SCM and GLCM descriptors. As noted previously for the detection of vineyards, the computational cost of the maximum likelihood classification with waveledbased multivariate models is as expected significantly higher than the one observed for the k-NN classifier with GLCM descriptors.

C. Discussion

In both applications detailed above, at least one of the tested wavelet-based multivariate models revealed higher classification accuracies than the use of GLCM descriptors, a standard approach for textural analysis in VHR optical remote sensing data. It was however observed that the models achieving the highest accuracy can be different depending on the context of application. The SCM and GGC models were indeed more accurate to detect vineyards while the SIRVg0 and GGC models displayed better performances in the case of oyster fields mapping. In a recent work [START_REF] Regniers | Wavelet-based texture features for the classification of age classes in a maritime pine forest[END_REF], a similar study was carried on the classification of mono-specific forest stands and showed that the SIRVgauss model was the best choice for this particular application. Each context of application has indeed its own specificities in terms of size, frequency and anisotropy of the textures to differentiate. The best model to be used can thereby be different depending on these specificities. A CBIR analysis carried on the learning database has proven to be a good way to identify the most adequate model prior to the classification.

Overall, this confirms the interest of the use of multivariate models to characterize textures and the potentially richer textural information captured through this modeling approach in the wavelet domain. The proposed supervised classification strategy relying on the use of a learning database composed of texture-wise homogeneous patches or ROI and a presegmentation of the image to classify also revealed to be an appropriate procedure and achieved satisfactory classification results in both applications. This strategy was applied by using a DWT, this choice being motivated by the ability of the DWT to produce decorrelated subbands. Nonetheless, any other types of wavelet decompositions (stationary wavelet transform, complex wavelet transform, steerable pyramids, etc.) could be used. Besides, one major quality of the proposed strategy is that only a few parameters need to be set. The only choices to be made concern the type of wavelet decomposition and the number of scales of this decomposition, as well as the eventual parameters related to the segmentation algorithm.

Once these options are set, the rest of the procedure is straightforward and can be easily adapted to any kind of texture-based classification problems.

It is also worth noting that texture features derived from DWT wavelet subbands are specific to the orientation of each subband and features are not in essence rotation invariant. Therefore, two identical anisotropic textures only differing by their main orientation could be considered as different or dissimilar when using such features. This problem comes from the similarity measure used between models as the distance is measured between pairs of orientation specific subbands. In the applications exposed above, the lack of rotation invariance was compensated by using a fairly exhaustive learning database in which all possible orientations of the textures are represented. In other cases, this exhaustivity can be more difficult to reach. In [START_REF] Regniers | Wavelet-based texture features for the classification of age classes in a maritime pine forest[END_REF], a solution was proposed by applying a rotation to each region of the learning database and to each region to classify in order to force the main orientation to 0°before proceeding with the feature extraction. This solution proved to be efficient with a significant increase in classification accuracy.

Finally, we employed here a Mean-Shift algorithm to perform the pre-segmentation to identify homogeneous regions in the image to classify. As it was observed in the two presented applications, this algorithm has some drawbacks, especially in narrow areas where regions of mixed textures can be created. A pre-treatment filtering process could be used prior to the segmentation to avoid these mixed texture regions by applying for example anisotropic diffusion filters [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] or by using morphological neighborhoods [START_REF] Fauvel | A spatial-spectral kernel-based approach for the classification of remote sensing images[END_REF]. Besides, any other types of segmentation algorithm could be used to produce a pre-partition of the image.

V. CONCLUSION

Although wavelet-based probabilistic models gained in popularity in the image processing and computer vision communities in the recent years, they have not yet been extensively explored for applications in VHR optical remote sensing data. Here, we gathered evidence that this approach is also perfectly suitable for the supervised classification of land covers in such VHR optical data. Moreover, we proposed a complete strategy to apply wavelet-based multivariate models in a supervised classification procedure of textures in VHR panchromatic data. This strategy relies on a learning texture database made of texture patches or ROI and on a pre-partition of the image to classify. Texture features are extracted from the learning database and from the regions of the pre-partition by using multivariate models (SCM, SIRVgauss, SIRVg0, GGC) to represent the distribution of observed local spatial dependencies in wavelet subbands in a multi-scale and multi-orientation framework. A CBIR analysis carried on the learning database is first conducted to identify the most efficient models to retrieve textures in the context of application. A classifier based on a similarity measure or a likelihood criterion is next used to produce classification results with the most performant models. The applicability of this strategy was tested in two distinct contexts. In both applications, the use of the proposed strategy has enabled to achieve satisfactory classification results with at least one of the tested multivariate models displaying higher classification accuracies than the standard texture analysis using GLCM descriptors. These results confirm the pertinence of using multivariate modeling in the wavelet domain to capture potentially complex texture patterns and improve the classification of VHR optical remote sensing data. Moreover, the classification is almost straightforward with only a few parameters to be set. In terms of perspective, this strategy and the use of multivariate models could also be tested in a context of multi-sensor data fusion by combining the modeling of dependencies between color bands in mutispectral data or polarimetric bands in SAR data with the modeling of spatial dependencies in panchromatic image.

Fig. 1 :

 1 Fig. 1: Modeling of discrete wavelet subbands where k i are the random vectors of observation, η are the parameters of the model P , N sc and N or are respectively the number of scales and orientations of the decomposition.

Fig. 3 :

 3 Fig. 3: Examples of texture patches used in a learning database for vineyards classification.

Fig. 4 :

 4 Fig. 4: Recall/Precision curves for GLCM and four multivariate models obtained with vineyard texture databases created from four panchromatic Pléiades data.

Fig. 5 :

 5 Fig. 5: Classification results on an extract of the 2012-08-22 Pléiades data acquired in Pessac-Léognan appellation. Green contours in Fig. 5b correspond to the frontiers of the regions obtained with the Mean-Shift algorithm.

Fig. 6 :

 6 Fig. 6: Examples of texture ROI used in the learning database for oyster fields classification.

Fig. 7 :

 7 Fig. 7: Recall/Precision curve for GLCM and four multivariate models obtained with the texture database created from 2013-04-25 panchromatic Pléiades data.

Fig. 8 :

 8 Fig. 8: Classification results on an extract of the 2013-04-25 Pléiades data acquired on the Arcachon Bay. Green contours in Fig. 8b correspond to the frontiers of the regions obtained with the Mean-Shift algorithm.

TABLE I :

 I Accuracy indices and computational time for vineyards classification

Table II .

 II These results are discussed here below.

	Method	OA	kappa	PA	UA	Time 1
	GLCM -kNN 86.28	0.66	83.13	69.2	47 sec.
	SCM -ML	89.39	0.61	82.88 61.23	98 min.
	SIRVg0 -ML	86.3	0.7	87.72 71.44	104 min.
	GGC -ML	88.01	0.7	84.27 72.87	109 min.
	OA = Overall Accuracy -PA = Producer's Accuracy -UA = User's
			Accuracy

TABLE II :

 II Accuracy evaluation for oyster racks classification

ACKNOWLEDGMENT

The authors would like to thank the Recette Thématique Utilisateur/ORFEO and ISIS program of the CNES for providing the Pléiades images used in this study.

This work was supported by the Aquitaine Regional Council and CNES (Centre National d' Études Spatiales) through the TOSCA project STELLA THR. O.


Yusuf Bhujwalla
email: yusuf-michael.bhujwalla@univ-lorraine.fr

Vincent Laurain
email: vincent.laurain@univ-lorraine.fr

Marion Gilson
email: marion.gilson@univ-lorraine.fr

The Impact of Smoothness on Model Class Selection in Nonlinear System Identification: An Application of Derivatives in the RKHS

Keywords: Nonlinear System Identification, Nonparametric Modeling, RKHS, Regularization, Derivatives, Gradient Regularization

The impact of smoothness on model class selection in nonlinear system identification: An application of derivatives in the RKHS

I. INTRODUCTION

Many real systems are nonlinear, or exhibit some nonlinear behavior. As such, nonlinear identification represents a wellstudied and important area of system identification [START_REF] Billings | Identification of non-linear rational systems using a prediction-error estimation algorithm[END_REF], [START_REF] Nelles | Nonlinear System Identification: From Classical Approaches to Neural Networks and Fuzzy Models[END_REF]. However, attempting to estimate nonlinear models can be very challenging, posing many problems of which the user must be aware [START_REF] Sjöberg | Nonlinear black-box modeling in system identification: A unified overview[END_REF].

One problem is that nonlinear models represent a vast range of functions. Defining a model class that encompasses the true behavior of a system can be extremely difficult. In the absence of sufficient prior knowledge about the system, how can any useful statements about these nonlinearities be made?

A way of dealing with this is to consider a nonparametric representation of the system. Here, the number of parameters describing a model depends on the amount of data available, offering a very flexible way of describing nonlinear systems. However, this flexible over-parameterized approach has certain drawbacks. For a model of this form, its estimated parameters will typically have unacceptably high variance, and be a very poor predictor of future outputs. This is commonly dealt with using regularization.

Regularization is an established statistical technique [START_REF] Chen | On the estimation of transfer functions, regularizations and gaussian processes -revisited[END_REF], [START_REF] Sjöberg | On the use of regularization in system identification[END_REF], [START_REF] Tikhonov | Solutions of Ill-Posed Problems[END_REF], allowing the incorporation of some prior intuition of how the true system may behave into the model. It is typically used as a way of controlling the bias-variance trade-off in a model, in variable selection, or as a way of enforcing some desired properties upon the estimator (such as smoothness [START_REF] Duijkers | Shrinking complexity of scheduling dependencies in ls-svm based lpv system identification[END_REF] or sparsity [START_REF] Rosasco | Nonparametric sparsity and regularization[END_REF]).

An elegant way of constructing a nonparametric identification problem is to use the theory of Reproducing Kernel Hilbert Spaces (RKHS). It should be mentioned that there is an equivalence between RKHS methods and several other branches of kernel methods, notably the Gaussian Process and Least-Squares Support Vector Machine approaches in machine learning. This framework is applicable to many different model classes, and allows their implicit definition through an associated kernel function.

Whilst the regularization and its associated hyperparameters are often discussed, the choice of a suitable kernel is frequently overlooked. As the kernel defines the RKHS, it also implicitly defines the model class -and hence the capacity of the model to capture the true nonlinearity and the smoothness of the true function.

In this paper, we will investigate if, by constraining the smoothness of the estimated function through the regularization, it is possible to integrate the kernel selection into the optimization criterion, instead of relying on an a priori choice of kernel. This paper will be structured as follows. In section II, the identification problem will be formulated. Section III discusses the smoothness-flexibility trade-off, and the impact this has on the identification problem. In section IV, a potential solution to this problem is proposed, along with some results regarding derivatives in the RKHS. Finally, in section V, a simulation example will be given, providing an indication of the performance of the proposed scheme and hopefully illustrating its potential advantages.

II. PROBLEM DESCRIPTION

A. The Data-Generating System

We will now introduce the identification setup. Assuming we have N observations of a data-generating system S o :

D N = {(u 1 , y 1 ), (u 2 , y 2 ), • • • , (u N , y N )} . (1)
In System Identification, often a variable x k ∈ R na+n b +1 is introduced. Typically, x k is the regressor vector, composed of the past and present inputs and outputs of the system, u and y:

x k = [y k-1 • • • y k-na u k • • • u k-n b ] (2)
where n a and n b are the orders of the output and input respectively. Using these definitions, we will attempt to reconstruct a general nonlinear model of the unknown true system S o describing its behavior:

y k = f o (x k ) + e o,k (3)
Here, f o : R na+n b +1 → R is an unknown nonlinear function and e o,k is an additive noise term at the output at each state k. To simplify the analysis, here S o is assumed to be a single-input single-output (SISO) system, but the extension to the multivariate case is straightforward.

It will also be assumed that e o is white Gaussian noise (see [START_REF] Laurain | An instrumental least squares support vector machine for nonlinear systems identification[END_REF] for a discussion of the effect of colored noise in identification).

B. Hilbert Spaces

The function f (•) is stated to be part of the Hilbert space H (f ∈ H : X → R). A Hilbert space represents all possible realizations of some particular class of functions, for example all functions of continuity degree C k . Moreover, a Hilbert space is a vector space such that any function f ∈ H must have a nonnegative norm, f H > 0 (for f = 0). f must be also equipped with an inner-product in H. For example, for real-valued functions f, g ∈ L 2 , the inner-product would have the familiar form:

f (x), g(x) = X f (x)g(x) dx. (4)
The properties of Hilbert spaces have been explored in great detail throughout the literature. For the interested reader, it is recommended to refer to one of the many texts discussing the subject, such as [START_REF] Zhou | Robust and Optimal Control[END_REF].

C. Reproducing Kernels

An extremely useful property of Hilbert Spaces is their equivalence with an associated kernel function [START_REF] Aronszajn | Theory of reproducing kernels[END_REF]. This equivalence allows us the simple definition of a kernel, instead of fully defining the associated vector space. More formally, if a Hilbert space H is a Reproducing Kernel Hilbert Space (RKHS), it will have a unique kernel, K : X × X → R, spanning the space H.

To simplify the analysis, it will be assumed that x k = u k , and therefore our observations of the system are now {x i } N i=1 ∈ R N . Importantly, any function in H can be represented as a infinite weighted linear sum of this kernel evaluated over the space H, as:

f (•) = f (x), k x (•) H = ∞ i=1 α i k xi (•) (5)
where {α i } ∞ i=1 ∈ R ∞ are the parameters of the model. This relationship is known as the Reproducing Property and k xi is known as the kernel slice.

One of the most useful properties of a kernel in the RKHS is that:

k xi , k xj = k xi (x j ) = k xj (x i ) = K(x i , x j ) (6)
that is, the dot-product of a reproducing kernel with itself is itself the kernel. This result is commonly known as the kernel trick, and allows us to write the inner-product as a tractable function which implicitly defines a higher (or even infinite) dimensional space.

For K to be a valid kernel in the RKHS, it must be a Mercer kernel [START_REF] Schölkopf | A generalized representer theorem[END_REF].

D. The Representer Theorem

In reality any estimation problem will deal only with finite data. But since f is represented by infinitely many parameters, a regularized cost function is commonly used in the form:

V(e, f ) = c((x 1 , y 1 , f (x 1 )), • • • , (x N , y N , f (x N ))) + g( f H ). (7)
For example, c(•) is often a loss-function in L 2 , such that

c (x, y, f (x))) = N k=1 (y k -f (x k )) 2 (8)
where

x = [x 1 • • • x N ] and y = [y 1 • • • y N ] . g(•)
is an additional global constraint on the function, independent of the observations. This term is the regularization term. Often, g(•) is taken as a constant:

g ( f H ) = λ f H (9)
Here λ is the regularization hyperparameter -which can be considered as controlling the bias-variance trade-off. If f is in H, (7) permits a truncated form of the expression given in [START_REF] Kimeldorf | Some results on tchebycheffian spline functions[END_REF]. This means the infinite expansion given in [START_REF] Kimeldorf | Some results on tchebycheffian spline functions[END_REF] reduces to a finite summation around the observations alone. This is The Representer Theorem, for which a proof can be found in [START_REF] Schölkopf | A generalized representer theorem[END_REF]. For V(e, f ) given above, f can be written as

f (•) = N i=1 α i k xi (•), i = 1, 2, . . . , N α i ∈ R (10)
Using this expression, the norm f H , can now be defined as

f H = f, f H = N i=1 N j=1 α i α j K(x i , x j ) = α Kα ( 11
)
where

K is the N × N Gram Matrix, {K} i,j = K(x i , x j ) and α = [α 1 • • • α N ] is the parameter vector.
Taking ( 8) and ( 9) yields the following cost-function:

V(e, f ) = y -f (x) 2 2 + λ f H . (12)
A closed-form solution can now be obtained for α:

⇒ α = (K + λI) -1 y. (13)
This solution will likely be familiar to the reader, as it is widely used in inverse problems. Note that this solution will depend on the hyperparameter λ and the choice of kernel K.

Typically, determination of hyperparameters is a non-trivial problem, requiring the implementation of methods such as cross-validation or marginal likelihood [START_REF] Pillonetto | Kernel methods in system identification, machine learning and function estimation: A survey[END_REF], [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF].

III. MODEL CLASS SELECTION A. The Kernel Selection Dilemma

Whilst the choice of λ is often discussed in the literature, the choice of the kernel K is usually left somewhat open. But in fact, this is a very important part of the identification process. Selecting a kernel usually requires some statement about the nature of the system in question to be made.

A range of kernel functions are available to the user, for example step, linear, polynomial and spline functions. A popular choice of kernel for continuous functions, which we will consider here, is the Gaussian RBF kernel:

k xi (x) = exp - x -x i 2 σ 2 . ( 14
)
which represents the C ∞ class of functions. Henceforth, the discussion will be continued solely considering RBF kernels.

The results developed will nonetheless be valid independent of the choice of kernel.

In the RBF case, the hyperparameter σ defines the width of K, and hence the smoothness of the estimated function. Increasing σ can be considered as smoothing f by progressively filtering out higher-frequencies. As σ characterizes K, and K defines H, σ also influences H. Hence, σ is somehow related to the model class selection. However, the choice of σ will also be influenced by the noise. We will now introduce an example to illustrate this paradox.

B. A Simple Example

A one-dimensional static switching function is corrupted by white noise at the output. A signal-to-noise ratio of SNR = 10dB is used in this example, defined by the ratio of the noise-free output to the noise level: SNR = 20 log σ ẙ/σ e . No other information about the function, such as the switching points, is provided.

The unknown system S o is described by the following equation:

y k = f o (x k ) + e o,k , where x k ∼ U (-1, 1
) is a uniformly distributed input-signal and e o,k ∼ N (0, σ 2 N ) is white Gaussian noise. The function f o is a non-zero mean, non-smooth function:

f o (x) = 15,
x ≤ -0.5, or 0 < x ≤ 0.5 -5, 0.5 < x ≤ 0, or 0.5 < x

In Figure 1, each image shows the function reconstructed for different values of σ and λ. As mentioned, an RBF kernel was used, as given in [START_REF] Sjöberg | On the use of regularization in system identification[END_REF]. Clearly, the discontinuous f o ∈ C -1 lies outside the Hilbert space of the RBF kernel -as it would require an infinite range of frequencies to fully reconstruct the function. Using such a function allows us to examine the effect of σ and λ on the limits of the available model class.

The kernel width needs to be sufficiently small not to make any assumptions on the smoothness of f , allowing a large model class. In the top-left figure, it can be seen that choosing such a σ allows the dynamics of the true function to be captured by the estimated function -but it is unacceptably sensitive to noise.

In the bottom-left figure, we see that increasing σ does improve the variance characteristics of the estimate, but at the cost of placing a strong assumption on the smoothness of the function. The effect of reducing the size of the available model class means we have failed to capture the dynamics of the system.

From [START_REF] Pillonetto | Kernel methods in system identification, machine learning and function estimation: A survey[END_REF], it is possible to introduce a regularization term to improve the characteristics of the estimator. But, as can be seen in the top-right and bottom-right figures, even when an excessively high λ is used to reduce the variance, such that an undesired bias is introduced, the function is still not smooth.

Optimal values for the hyperparameters can be found, but this does not change that, in practice, we are forced to choose between the flexibility of our model and the smoothness.

Hence, when selecting a kernel, how can we choose the widest class possible in order to avoid making assumptions on f o , whilst still having a reliable estimator, capable of reproducing any given data and predicting future outputs?

In the following section, we will propose a method to directly enforce smoothness independently of the model class selection by constraining the gradient of the estimated function.

IV. SMOOTHNESS IN THE RKHS

The smoothness of a function f is inextricably linked to the continuity of the space C k within which it resides. Therefore, a natural way of enforcing smoothness on a function would be to place constraints on its derivatives:

df (x) dx , d 2 f (x) dx 2 , • • • , d k f (x) dx k (15)
There are several ways this could be achieved. Before proposing a solution in section IV-C, we will briefly review several existing methods from the literature.

A. Sobolev Spaces

One relatively well-known method of enforcing smoothness involves the consideration of Sobolev Spaces. In a Sobolev space, S k,p of continuity C k and norm p, a function within the space is quantified in terms of its derivative. For p = 2, this Sobolev Space becomes a Hilbert Space, H k , and we can define:

f H k = k i=0 X d i f (x) dx i 2 dx (16)
Using this definition of the norm will naturally impose continuity constraints on the function up to order k. Under certain conditions, a minimal representation of f in this case can be found based on the spline kernel, as discussed in [START_REF] Wahba | Spline models for observational data[END_REF]. Though this is a very effective solution, we would like to develop a method valid for a range of kernels, rather than enforcing the choice of kernel.

B. Identification with Observations on the Derivative

Another method discussed in the literature is based on minimization against observations of the derivative [START_REF] Rosasco | A regularization approach to nonlinear variable selection[END_REF]. In certain cases it may not be feasible to obtain measurements directly for the output, or it may be of interest to constrain the error against the derivatives. Then, a cost function of the form below may be used:

V obvs (f ) = y -f (x) 2 2 + γ 1 dy dx - df (x) dx 2 2 + • • • γ m d m y dx m - d m f (x) dx m 2 2 + λ f H (17)
where here we have considered an MSE loss. Again, this is capable of enforcing smoothness, but it introduces m additional hyperparameters (denoted here by γ) increasing the complexity of the problem.

C. Regularization using Derivatives

Here, a scheme independent of kernel choice is proposed, with the advantages of neither requiring any additional hyperparameters nor any significant additional computational load (compared with the solution given in ( 13)). By replacing the regularization term given in [START_REF] Ljung | System Identification, theory for the user[END_REF] with one optimizing against the derivative of the function, we have:

V ∇ (e, f ) = y -f (x) 2 2 + λ Df H ( 18
)
where

D m f = d m f (x)
dx m is the m th order differential operator. Although not identical, there is a clear parallel with the spline-smoothing problem mentioned in section IV-A. Now λ controls the smoothness to be put on the function, meaning σ is no longer critical in reducing the sensitivity of the model to noise. By explicitly forcing the regularization to act in this way, a greater flexibility in the model class may be permitted.

Note that as λ tends to ∞, for a bounded input this enforces a bounded solution, as f H → 0 now becomes Df H → 0, which implies that ∀x ∈ R:

⇒ lim λ→∞ f (x) = c, where 0 ≤ c < ∞ ∈ R (19)

D. The Representer Theorem for Derivative Regularization

As discussed in section II, for the cost-function given in [START_REF] Ljung | System Identification, theory for the user[END_REF], a valid representer is f (•) = N i=1 α i k xi (•). However, this form is not necessarily valid for the case of [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF].

For the proof of the generalized representer theorem given in [START_REF] Schölkopf | A generalized representer theorem[END_REF] to hold, it is required that a relationship can be defined such that

D m f H = g ( f H ) ( 20
)
where g(•) is a strictly monotically increasing function on the norm of f . Whilst a relationship can be defined, it will not in general be a monotonic function of f H . In the case of the Gaussian RBF kernel adopted in section III-B, (20) does not hold. Hence, a representer of the form stated above will be suboptimal for [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF]. This means that, unlike in the case of ( 7), the reproducing property of ( 5) cannot strictly be truncated to a finite expression evaluated solely over the observed data.

Nonetheless, we will proceed using this form, whilst acknowledging that a better-performing solution may exist. Whilst using a suboptimal representer is not ideal, it will allow us to preserve the mathematical simplicity and computational efficiency of the solution given in [START_REF] Schölkopf | A generalized representer theorem[END_REF]. Furthermore, in addition to being able to control the smoothness of the estimated function, for certain choices of kernels the statistical properties of (7) are loosely preserved:

Lemma 1 (The boundedness of f (x) for ( 18)). For a costfunction of the form given in [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF], let the representer of the estimated function be f

(x) = N i=1 α i k xi (x). If f (x) : X → R ∈ H is bounded such that f H < ∞, and the kernel K : X ×X → R spanning H decays such that k xi (x) -→ 0 ||x-xi||→∞ , f (x) is bounded as: lim λ→∞ f (x) = 0, ∀x ∈ R. (21)
Proof:

Given that f H < ∞ and k xi (x) → 0 as x -x i → ∞: ⇒ lim x-xi →∞ f (x) = 0 (22)
where [START_REF] Zhou | Robust and Optimal Control[END_REF], that for [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF]:

{x i } N i=1 ∈ X N , ∀i = 1 • • • N . Now recall from
lim λ→∞ f (x) = c, ∀x ∈ R (23)
where 0 ≤ c < ∞ ∈ R.

Observe that (23) applies across R. Hence if ∃x ∈ R as in ( 22), c must be 0. Therefore:

lim λ→∞ f (x) = 0, ∀x ∈ R. (24)
Hence, under the above-stated conditions, f (x) is bounded by the regularization term of (18) as in [START_REF] Ljung | System Identification, theory for the user[END_REF]. This suboptimal representer thus allows the implicit regularization of f H in [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF], in addition to enforcing smoothness on f (x).

E. A Closed-Form Solution

From [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF], let f ∈ H k , where k denotes the order of continuity of the space H k . Now, for x ∈ X where again we consider X = R, D m f (x) ∈ H k provided m < k. This allows a reproducing property to be defined for derivatives of functions, analogous to that stated in [START_REF] Kimeldorf | Some results on tchebycheffian spline functions[END_REF]:

D m f (•) = f (x), D [0 m] k x (•) H = ∞ i=1 α i D [0 m] k xi (•) (25)
where the operator D [ι κ] is defined as differentiation with respect to the first and second variables of the kernel function as:

D [ι κ] k xi (x) = d ι+κ k xi (x) dx κ dx ι i . (26)
Now we can proceed to define D m f H . From [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF], 18), m = 1, giving:

D [0 m] K, D [0 m] K H = D [m m] K. In (
Df H = Df, Df H = N i=1 N j=1 α i α j d 2 k xi (x j ) dx j dx i = α D [1 1] Kα (27)
where

D [1 1] K is the N × N matrix of derivatives {D [1 1] K} i,j = d 2 k xi (x j ) dx j dx i .
In the RBF case, the derivative is given by:

d 2 k xi (x) dx dx i = 2 σ 2 1 - 2 σ 2 (x -x i ) 2 exp - x -x i 2 σ 2 .
(28) From ( 27), the derivation of a closed-form solution for α ∇ , the optimal parameters of V ∇ , is trivial:

⇒ α ∇ = K K + λD [1 1] K -1 K y (29)
Note that this solution has a similar computational complexity to the solution given in [START_REF] Schölkopf | A generalized representer theorem[END_REF].

F. A Simple Example

We now demonstrate the effect of regularizing over the gradient on the example from section III-B.

Figure 2 shows the estimation function for a small kernel (σ = 0.01) over a range of λ values with comparable variances in each case. The effects of the regularization term in the two cases are drastically different. On the left-hand side, the effect of λ on the bias-variance trade-off is clearly visible; the variance is progressively reduced, at the cost of an increasing bias in the model. But, the regularization does not act on the smoothness of the function.

On the right-hand side, the variance is reduced by smoothing the function. Eventually the estimated function becomes biased and excessively smooth. In practice, this simply means that λ must be tuned properly.

V. SIMULATION EXAMPLE

A Monte-Carlo study based on a simulation example will now be used to compare the two algorithms discussed in this paper -the f H -regularized and Df H -regularized approaches.

A. The Data-Generating System

Two 1-D nonlinear functions are corrupted with white Gaussian noise at the output (SNR = 10dB). Both systems were excited using N = 1000 uniformly distributed datapoints, x k ∼ U (-1, 1). The first system S 1 o , a smooth function, is given by the following equation:

f 1 o (x) = 10((1 -x) -2sinc(4x) 2 , (30)
and the second system S 2 o , a non-smooth function, by:

f 2 o (x) = 10(1 -x), |x| > 0.2 10 [(1 -x) + 10(|x| -0.2)] , |x| ≤ 0.2
Two approaches to the hyperparameterization were used, with 100 Monte-Carlo trials run in each case. For both approaches, the kernel function was chosen to be the Gaussian RBF kernel defined in section III. 1. For f H -reg: λ and σ were optimized using crossvalidation on a separate noisy validation dataset.

2. For Df H -reg: the kernel width was fixed a priori to σ = 0.01 such that the smallest value capable of enforcing smoothness is chosen (σ ∼ max(x i+1 -x i ), i = 1 . . . N -1). λ was determined by cross-validation.

B. Estimation Results

Figure 3 shows the results of the Monte Carlo simulation. The average estimated function and the variance of the output across the input space are plotted for both methods. As an illustration, a sample kernel is given in each case. The central column shows the true nonlinearity and a sample of the noisy estimation data for both functions.

Table 1 shows the hyperparameter values obtained by cross-validation in each case. In addition, the variance of the estimator averaged of the input space and the mean best fit rate against the noise-free function (BFR) over the trials are given (where BFR = 100 1 -ŷ -y 2 y -ȳ 2 ).

TABLE I SUMMARIZED RESULTS

From Table 1 it can be seen that both approaches estimate the function quite succesfully. However the f H -Reg scheme outperforms the Df H -Reg scheme. Despite this, analysis of Figure 3 shows several advantages of the Df H -Reg approach.

Firstly, using a small kernel, we can still estimate the smooth S 1 o -even with large levels of noise. In the f H -Reg approach, a much broader kernel is required. But, as for the Df H -Reg a small kernel can be chosen, we are never forced to select an overly smooth model class. As such, the structural difference between the smooth S 1 o and nonsmooth S 2 o is clearly apparent. For f H -Reg, it is not. Furthermore, the results of Df H -Reg were achieved without having to optimize over the kernel hyperparameter. For X = R, this is already advantageous. As the dimensionality of the input increases, this becomes increasing attractive.

VI. CONCLUSIONS

In this paper, a derivative-based regularization method has been proposed, attempting to divorce the dependency of the smoothness of the estimated function from the selection of the underlying model class. It has been shown that using this scheme, both smooth and non-smooth functions can be estimated by using a small kernel width and shifting the hyperparameterization problem solely onto λ.

With λ appearing linearly in the parameters, this may simplify the hyperparameter optimization or perhaps permit the introduction of a varying smoothness constraint across the input space. Consideration of the multi-dimensional case and further investigation of the sub-optimality of the representer are intended as future works.

Fig. 1 .

 1 Fig. 1. Estimation of 1D switching signal for different hyperparameter values. The mean value of each function over 1000 Monte-Carlo trials is plotted against the true function, with ± the variance marked on either side.

Fig. 2 .

 2 Fig. 2. Monte-Carlo Results: Mean Estimates ± variance for a small kernel (σ = 0.01), with λ values given on the y-axis of each plot

Fig. 3 .

 3 Fig. 3. Mean estimates of the functions, ± variance, and the kernel size used in each case for the smooth f 1 (top) and non-smooth f 2 (bottom); f H -Reg (left) and Df H -Reg (right)