B Chalmond
email: bernard.chalmond@cmla.ens-cachan.fr

S Girard
email: girard@stat.math.univ-montp2.fr

Non linear modelling of scattered multivariate data and its application to shape change

Keywords: Data analysis, Example-based analysis and synthesis, Visual learning, Face representation, Principal components analysis, Nonlinear PCA models, Dimensionality reduction, Multidimensionnal scaling, Projection pursuit, Eingenfeatures

We are given a set of points in a space of high dimension. For instance, this set may represent many visual appearances of an object, a face or a hand. We address the problem of approximating this set by a manifold in order to have a compact representation of the object appearance. When the scattering of this set is approximately an ellipsoid, then the problem has a well-known solution given by Principal Components Analysis (PCA). However, in some situations like object displacement learning or face learning this linear technique may be ill-adapted and nonlinear approximation has to be introduced. The method we propose can be seen as a Non Linear PCA (NLPCA), the main difficulty being that the data are not ordered. We propose an index which favours the choice of axes preserving the neighborhood of the nearest neighbours. These axes determine an order for visiting all the points when smoothing. Finally a new criterion, called "generalization error", is introduced to determine the smoothing rate, that is the knot number of the spline fitting. Experimental results conclude this paper: the method is tested on artificial data and on two data bases used in visual learning.

Introduction

In many situations in pattern recognition, machine intelligence, or artificial vision, it is necessary to approximate multivariate data by a parametric model in order to be able to handle more easily the information contained in these data. In this paper, the data set X = {x j , 1 ≤ j ≤ N } belongs to a high dimensional space : x j = (x (i) j ) 1≤i≤n ∈ IR n . The problem is then to compress this data set to a low-dimensional manifold. A technique now commonly used for dimensionality reduction in computer vision is Principal Components Analysis (PCA) which yields a linear representation. Since this technique summarizes the data by the mean and the standard deviation (the covariance matrix), the linear representation is accurate only if the data distribution is gaussian, i.e., if the cloud X is a n-D ellipsoid (Fig. 2(a)). Although this linear model is effective in a very broad range of applications, there are however situations where the PCA breaks down. It happens when means and standard deviations do not accuratly reflect the data distribution, in other words, when the cloud is not an ellipsoid (Fig. 2(c)). This indicates the need for nonlinear representations.

In this paper a novel technique for designing and fitting nonlinear models is proposed.

From the geometric point of view, we try to approximate X by a d-dimensional manifold with d < n [START_REF] Milnor | Topology from the differentiable point of view[END_REF]. This manifold will be given by an implicit equation G(θ, x) = 0, where θ is a vector of parameters and G(θ, .) is a continuously differentiable function from IR n to IR n . A model is a particular analytical expression of G, the most simple being the linear representation as it is given by Principal Components Analysis. Our problem, in a nutshell, is a smoothing problem. However in our case, it is not a classical smoothing problem since the points x j are not ordered. Fitting a manifold is a much more difficult problem than function fitting. A second difficulty adds to the former. The number of points x j in the data set is usually moderate with respect to the dimension n. This means that the space IR n is nearly empty, what it is known as the "curse of dimensionality" [START_REF] Huber | Projection Pursuit[END_REF].

Application domains

Let us illustrate our purpose by examples in visual learning, although our technique is not limited to this domain. The data set X is composed of N instances of an object appearance. Here, "object" is a generic term which designates for example a rigid or flexible manufactured object, a face, a hand,... [START_REF] Murase | Visual learning and recognition of 3D objects from appearance[END_REF][START_REF] Girard | Building and training radiographic models for flexible object identification from incomplete data[END_REF][START_REF] Girard | Auto-associative models, nonlinear Principal component analysis, manifolds and projection pursuit[END_REF][START_REF] Kirby | Application of the Karuhnen-Loève procedure for the characterisation of human faces[END_REF][START_REF] Kervrann | Learning structure deformation modes of nonrigid objects in long image sequences[END_REF]. It is convenient to see X as a set of occurrences of a random vector X in IR n . Such an occurrence can be a curve [START_REF] Girard | Designing nonlinear models for flexible curves[END_REF], an image [START_REF] Kirby | Application of the Karuhnen-Loève procedure for the characterisation of human faces[END_REF], or a feature vector as Gabor coefficients for example [START_REF] Mckenna | Face tracking and pose representation[END_REF]. An approach for the visual learning problem consists in fitting a compact model of the object's appearance to X . The aim of this appearance-based representation is to capture the flexibility or the changes of the object in a compact manner such that every occurrence of X can be approximated by a point of the manifold defined by the model. Such a representation has wide applications: recognition and pose of 3D objects [START_REF] Murase | Visual learning and recognition of 3D objects from appearance[END_REF], face recognition [START_REF] Kirby | Application of the Karuhnen-Loève procedure for the characterisation of human faces[END_REF][START_REF] Hallinan | A low-dimensional representation of human faces for arbitrary lighting conditions[END_REF][START_REF] Sung | Example-based learning for view-based human face detection[END_REF], face tracking [START_REF] Mckenna | Face tracking and pose representation[END_REF], gesture recognition, [START_REF] Pavlovic | Visual interpretation for human-computer interaction: a review[END_REF], image retrieval [START_REF] Swets | Using discriminant eigenfeatures for image retrieval[END_REF]. In visual learning, PCA is widely used for dimensionality reduction to enable efficient learning. However, PCA use is limited since its linear representation is inappropriate for modelling nonlinear effects such as bending or rotation of shape. This drawback has been highlighted by Murase and Nayar [START_REF] Murase | Visual learning and recognition of 3D objects from appearance[END_REF]. These authors have used the principal components of many views of a single rigid object to visualise the low dimensional manifold describing changes due to rotation and illumination conditions. The object's pose is then determined by the position of its projection on this manifold.

Let us say a few words about the application which has motivated our research. In [START_REF] Girard | Building and training radiographic models for flexible object identification from incomplete data[END_REF],

the authors address the problem of identifying the radiographic projection of an object (its appearance) from incomplete data extracted from a radiographic image. They assume that the unknow appearance is a particular instance of the projection of a flexible object. Their approach consists in learning a deformation model and its probability distribution able to represent and to simulate a great variety of appearances. This modelization is achieved using a training set of complete appearances. Then, given the incomplete data, the identification task consists in estimating the unknown appearance using the previous model whose probability distribution plays the role of a prior distribution in a Bayesian framework.

PCA and its limitations

Let us recall some important points of the PCA method [START_REF] Mardia | Multivariate Analysis[END_REF][START_REF] Girard | A nonlinear PCA based on manifold approximation[END_REF] which will be usefull for the forthcoming nonlinear model. The PCA model is linear. When the data are centered, its expression is

G(a, x) def = x -F (a, x) (1)
with

F (a, x) = d k=1 a k (a T k x) ,
where a k ∈ IR n are unit orthogonal vectors which span the principal axes. Then, G(a, x) = 0 defines a d-dimensional linear subspace in IR n . As the x j are observations of a non-degenerated random variable, it is clear that they cannot verify the deterministic equation G(a, x j ) = 0 all together : G(a, x j ) = 0 and we define

r j (d) def = G(a, x j ) . (2)
For a d-dimensional model, r j (d) is interpreted as the residual approximation error of x j by F (a, x j ). The model ( 1) is said to be auto-associative since x j is approximated by a function of {x j }. The model parameters a = (a 1 , ..., a d ) in (1) are estimated by minimizing the empirical mean square N j=1 r j (d) 2 . Since these estimations are related to the x j orthogonal projections on the desired subspace, the r j are the Euclidian distances to the subspace.

Let us consider an artificial data set X composed of translated curves (Fig. 1(a)). This set corresponds to a cloud in IR n located on a 1-dimensional manifold. We can imagine what this manifold looks like, by projecting X on the first principal plane computed by PCA (Fig. 1(d)). This kind of figure is known in the statistical litterature as "the horseshoe effect" [START_REF] Mardia | Multivariate Analysis[END_REF]. We get a similar phenomenon when X is composed of faces in rotation [START_REF] Mckenna | Face tracking and pose representation[END_REF]. On the artificial training set, the estimated linear model leads to non-realistic curve simulations and does not allow to reproduce the variations due to a translation (Fig. 1(c)). Here, the simulated curves are obtained by drawing random points on the principal axes. This very simple academic example is only presented to illustrate our purpose and we must not mistake it for a registration problem [START_REF] Ramsay | Functional data analysis[END_REF]. Even for a variation as simple as translation, the linear model given by PCA is ill-adapted. More complex exemples are showned in Fig. 6 and Fig. 7. Recently, in the context of visual learning for object representation, nonlinear problems have arisen several times, but no general solution has been proposed [START_REF] Cootes | Active shape models, their training and application[END_REF][START_REF] Murase | Visual learning and recognition of 3D objects from appearance[END_REF].

Related work

The reconstruction of n -1 dimensional manifolds in IR n , like curves in IR 2 or surfaces in IR 3 , seems at first glance close to our problem. In [START_REF] Hoppe | Surface reconstruction from unorganized points[END_REF], a numerical algorithm is proposed. Its input is an unordered set X ⊂ IR 3 near an unknown manifold, and its output is a surface that approximates the manifold. Reconstruction problems of this kind occur in engineering applications like representation from range data. This algorithm makes heavy use of the fact that the data set is "continuously dense" in some sense, around the unknown manifold. This situation is not verified in our case. Furthermore, we are firstly interested in manifolds whose dimension d is much smaller than n, and secondly in representations based on an analytical model.

The PCA methods for analysing the variations of flexible curves, have grown in the last few years as a separate research topic called "Fonctional data analysis" [START_REF] Ramsay | Functional data analysis[END_REF]. The basic philosophy of functional data analysis is that one should think of observed data functions as single entities, rather than merely a sequence of individual observations. Nevertheless, the results of this research have not been generalized to flexible surfaces, i.e. to images.

Besides, PCA generalizations have been proposed in order to take into account nonlinear phenomena. PCA-like auto-associative methods have been studied from the neural networks point of view with perceptron networks [START_REF] Karuhnen | Generalization of Principal Component Analysis, optimization problems, and neural networks[END_REF], but in the end, these models remain linear. In the case n = 2, a truly nonlinear approach is proposed in [START_REF] Hastie | Principal Curves[END_REF]. It consists in searching a curve called "principal curve" which passes through the middle of the data set. This means that every point on the curve is the average of the observations projecting onto it. This yields a second approach to the n -1 dimensional manifofd reconstruction as introduced above. However, this technique is mainly dedicated to the dimension n = 2 and is non-parametric. Its extension to a greater dimension is a difficult task.

A specific technique close to PCA is Multimensional Scaling (MDS) [START_REF] Carroll | Multidimensionnal scaling[END_REF][START_REF] Kruskal | Multidimensional scaling[END_REF][START_REF] Shepard | Parametric representation of nonlinear data structures[END_REF]. MDS adresses the problem of constructing a configuration of

N points Y = {y j , 1 ≤ j ≤ N } in IR d satisfying y i -y j ≈ x i -x j .
This is achieved by minimizing a distorsion index whose expression is for instance [START_REF] Shepard | Parametric representation of nonlinear data structures[END_REF]:

I S (Y) = i<j ( y i -y j -x i -x j ) 2 i<j y i -y j 2 . ( 3
)
In this spirit, MDS technique has been applied to build self-organized neural networks which give a low-dimensional mapping of the manifod of a non linearly related data set. It yields an unfolding of the data manifold. We shall come back later to that kind of index.

In the second section of this paper, we define a Non Linear Principal Component Analysis (NLPCA) limited to one axis (d = 1) for sake of clarity. Our aim is to approximate X by a n-dimensional curve C. In other words, we want to determine a curve C which estimates an ideal curve C * using the sample X .

From this point of view, we have to deal with a curve reconstruction problem. Curve reconstruction has efficient solutions when the x j are ordered and a parametrization of the curve is given. So, as for MDS technique, our approach consists in searching a parametrization axis which yields an approximate ordering of the points x j . It follows that our methods works well when the curve to reconstruct is, roughly speaking, the graph of a function over an unknown axis in IR n . In IR 2 , this approach is less general than the principal curves one. In the third section, this is extended to several axis (d ≥ 1). Finally, in the fourth section, experimental results are described. 2 Nonlinear PCA with one axis

Introduction

Let us gradually introduce our algorithm in an intuitive fashion using Fig. 2. In Fig. 2(a), X is a non-isotrope Gaussian cloud and thus PCA is well adapted. Each point x j of X is approximated by its projection on the principal axis defined by the unit vector a 1 . This projection is:

F (a 1 , x j ) = a 1 (a T 1 x j ) .
In Fig. 2(b,c,d) the cloud X is no more Gaussian and thus its approximation by projection on the pincipal axis is inaccurate. In Fig. 2(b), we are naturally led to approximate X by a smooth curve with parameter b 1 . In that case, the x j approximation is

F (θ 1 , x j ) = S b 1 a T 1 x j with θ 1 = (a 1 , b 1 ) .
S is a function of u = a T 1 x. Let us note that u j = a T 1 x j is the coordinate of the x j projection on the axis a 1 . On this figure, the order of the u j in IR approximatively corresponds to the "topological order" of the x j in IR n . On the contrary, in Fig. 2(c) the principal axis a 1 does not allow to determinate the topological order. For instance, two non neighbor and distant points in IR n may have neighbor projections on the axis a 1 . It is a crucial problem. If we still attempt to approximate X by a smooth curve parametrized by u = a T 1 x, the representation would be inadequate for synthesis: outside the values u j , the function S b 1 (u) would not deliver realistic occurrences of X (Fig. 2(c)). This difficulty can be overcomed by choosing an axis which preserves the topological order after projection (Fig. 2(d)). For this, we shall define an index I(a, X ) measuring the preservation of the neighborhood structure of X . The selected axis will maximize this index :

a 1 = arg max a I(a, X ) (4)
Finally, the smooting function parameter is obtained as follows :

b 1 = arg min b 1 N N j=1 r j (1) 2 with r j (1) = x j -S b a T 1 x j . (5)
Our curve fitting algorithm consists of two stages: axis search (4) and smoothing [START_REF] Chalmond | Modelization for Image Analysis[END_REF]. The forthcoming sections develop in more details the successive steps of the algorithm.

Definition of I and S

In ( 4) and ( 5), we have to define the projection index I(a, X ) and the smoothing function S b .

Projection index

A first idea would be to use the index (3) which attempts to preserve all the distances. In our context, the important point is not the distance preservation but the preservation of the neighborhood structure of X : if two points are neighbors in IR n then their projections should be neighbors in IR and conversely, as it is illustrated in Fig. 2(d) but not in Fig. 2(c) for many points. To simplify this criterion, we decide not to preserve the complete neighborhood structure, but only the closest neighbors structure. The chosen index is naturally the number of points for which the constraints are verified:

I(a, X ) = N i=1 j =i Ψ[(x j closest to x i ) =⇒ (a T x j closest to a T x i )] , (6)
where Ψ is the indicator function. The analytical expression of I is obtained as follows. Let x φ(i) be the closest neighbor of x i . We want a T x φ(i) to be the closest neighbor of a T x i in IR.

Consequently, for every i, a has to verify N -1 inequalities:

I(a, X ) = N i=1 j =i Ψ a T x i -a T x φ(i) ≤ a T x i -a T x j = N i=1 j =i Ψ a T xiφ(i) ≤ a T xij ,
where xij = x i -x j . Finally, since a and -a define the same axis, we get:

I(a, X ) = N i=1 j =i Ψ a T xiφ(i) -xij ≥ 0 × Ψ a T xiφ(i) + xij ≤ 0 . (7)
It can be shown that this index has the following invariance properties: two invariance properties with respect to translation and scale indicate that this index belongs to the class III defined by Hubert [START_REF] Huber | Projection Pursuit[END_REF], which is well-adapted to Projection Pursuit algorithms.

The last property shows that the axis does not depend on the orientation of X (rotation and symmetry invariance).

Smoothing function

The

N } of the vector set X with respect to the parametrization {u 1 , ..., u N }, by a cubic regression spline

S (i) b (i) (u), u ∈ [u min , u max ].
To do that, we need to sort the set {u j } for each coordinate i. Let {u σ(j) } be the sorted set. The i-th coordinate smoothing is performed on the couples {(u σ(j) , x (i) σ(j) ), j = 1, ..., N }. Such a spline representation implies to choose discontinuity knots in the u σ(j) series. The number of knots tunes the smoothing rate. In general, this choice is quite difficult and will be discussed in Appendix.

Computation

Index maximization

To perform this optimization, we have to note that I(a, X ) is a finite piecewise constant function. Since the equation a T x = 0 implies that a belongs to the hyperplane orthogonal to x, the I(a, X ) expression [START_REF] Diaconis | On nonlinear functions of linear combinations[END_REF] shows N (N -1) hyperplanes noted H k :

H k = a ∈ IR n s.t. a T n k = 0 , 1 ≤ k ≤ N (N -1), n k ∈ xiφ(i) -xij , xiφ(i) + xij ; i = 1, ..., N ; j = i .
These N (N -1) hyperplanes determine a IR n partition in L regions R ℓ , (L ≤ 2 N (N -1) ), in which the index is constant:

I(a, X ) = L ℓ=1 I ℓ Ψ[a ∈ R ℓ ].
Although it is possible to compute I for every vector a by applying (7), we do not know how to extract a vector a from every region R ℓ . To overcome this difficulty, the idea is to construct an iterative algorithm A of type a(p + 1) = A[a(p)] which visits the regions from an initial solution a(0) (as a gradient algorithm would do it in the case of a regular function), a(p) and a(p + 1) being in distinct regions. The stochastic algorithms framework yields a general formalization to that optimization problem. To maximize I, the idea is to build a random walk upon the regions and try to increase the index value at each step. Each step is defined by a hyperplane H k obtained by random choice of the integers {k} and considered as a border of the new visited region. If a(p) is the solution at the previous step, the new candidate a(p + 1) is obtained by the orthogonal symmetry s k with respect to

H k : s k (a(p + 1)) = a(p) -2(a(p) T n k )n k . Algorithm A
1. Initialization: a(0) is chosen at random.

a(p) is known, a(p + 1) is computed as follows :

-k is randomly chosen between 1 and N (N -1).

-Comparison between the regions separated by H k :

a ′ ← s k (a(p)). If I(a ′ , X ) > I(a(p), X ) Then a(p + 1) ← a ′ Else a(p + 1) ← a(p).
3. Go back to 2, until convergence.

The resulting index sequence is decreasing, the algorithm converges quickly towards a local maximum. An iteration cost is equivalent to the index I(a ′ , X) computation cost, that is O(nN 2 ) (differences xij are computed as a preprocessing at the initialization step). Algorithm A ′ When a local maximum has been reached, all the regions accessible by the operator s k (dark grey regions in Fig. 3) have an index value lower than the maximum one (clear grey regions in Fig. 3). Yet, the index value on the other regions (white regions in Fig. 3) could be greater than the local maximum. As it is impossible to reach these regions directly from the region of local maximum, a solution for visiting them is to pass through a region of lower value (a dark grey region), the transition being governed by a probability. This is the basic idea of the simulated annealing that we briefly describe in our context [START_REF] Press | Numerical Recipes[END_REF]. Contrary to the algorithm A which imposes at each iteration p that ∆I =[I(a(p + 1), X ) -I(a(p), X )] > 0, the stochastic algorithm (algorithm A ′ ) authorizes ∆I to be negative: ∆I > T p log ξ where ξ is a uniform random number on ]0, 1[ and T p is a sequence decreasing towards 0 :

T p = T 0 λ p , with λ < 1.
The algoritm A ′ is the algorithm A in which the condition (If...Then) is modified by plugging ∆I > T p log ξ. As p grows T p log ξ → 0, and the probability to accept a region of lower index becomes zero. A rigorous writting of this algorithm can be derived from the Metropolis dynamic [START_REF] Chalmond | Modelization for Image Analysis[END_REF].

Spline smoothing

At this time, the principal axis is known through its estimation a 1 as obtained above. We have now to deal with a well-known smoothing problem since we use a spline regression model. Let us denote again u j = a T 1 x j the coordinate of the projection of x j . The goal here is to find b 1 which minimizes the approximation error

2 (b) = 1 N N j=1 r j (1) 2 = 1 N N j=1 x j -S b (u j ) 2 , = 1 N n i=1   N j=1 x (i) j -S (i) b (i) (u j ) 2   = n i=1 2 i (b (i) ). (8)
It appears that the approximation error 2 is expanded into n independent approximation errors 2 i . Let us note B the B-spline matrix corresponding to the chosen knots. Each of the n approximation errors is then written as

N j=1 x (i) σ(j) -S (i) b (i) (u σ(j) ) 2 = x (i) -Bb (i) 2 N , 1 ≤ i ≤ n , (9)
where . N is the IR N Euclidian norm. Let b

(i)
1 be the least-squares minimum, and ν the spline number of knots. Let us emphasize that b (i) 1 is computed as soon as ν is fixed, hence is denoted as b (i) 1 (ν). We use a second criterion to determine the best number of knots, as presented in Appendix.

Nonlinear PCA with several axes

In the spirit of the Projection Pursuit [START_REF] Huber | Projection Pursuit[END_REF], a second axis of projection can be computed on the residuals r j (1) which take the place of X in ( 4) and ( 5), leading to new residuals r j (2). This procedure can be again performed on the new residuals r j (2), and so on. NLPCA algorithm 1. Initialization: k ← 0 and r j (0) ← x j , j = 1, ..., N 2. Parameter estimation: θ k+1 = (a k+1 , b k+1 )

a k+1 = arg max a I(a, {r j (k)} N j=1 ) (10)
b k+1 = arg min b N j=1 r j (k) -S b (a T k+1 r j (k)) 2 (11)
3. Residual errors update:

r j (k + 1) ← r j (k) -S b k+1 (a T k+1 r j (k)) , j = 1, ..., N 4

Orthogonality conditions

We achieve the algorithm by adding some constraints that we have not given so far to simplify the presentation. Up to now, we have implicitly supposed that the residual errors r j (k) are orthogonal to the axis a k a T k r j (k) = 0 , ∀j .

It is straightforward to show that [START_REF] Friedman | Projection Pursuit Regression[END_REF] follows the natural condition a T k S b k (u) = u. Thanks to this condition the algorithm is very simple. After an approximation with respect to the axis a k , the residual errors r j (k) are located in its orthogonal subspace. So, the new axis as well as the new smoothing function are built in the a k , a k-1 , ..., a 1 orthogonal subspace, that is:

a T k a ℓ = 0 , ∀ 1 ≤ ℓ, k ≤ d (13)
a T ℓ S b k (u) = 0 , ∀ 1 ≤ ℓ < k ≤ d , u ∈ IR .
Let us emphasize that the orthogonality condition ( 13) gives an answer to the difficult problem of iterating the Projection Pursuit algorithm from a k to a k+1 , [START_REF] Friedman | Exploratory Projection Pursuit[END_REF].

Associated model

Now, the question is: what kind of representation is produced by our algorithm ? For d = 1, the answer is simple since X is approximated by a curve whose equation is [START_REF] Chalmond | Modelization for Image Analysis[END_REF]). In order to write this model for d > 1 in a concise form, let us denote a T k x = P a k (x), P a k being the function IR n → IR which defines the orthogonal projection of x j on the axis a k . With this notation, from the algorithm, we get:

G(θ 1 , x) = 0, G being the model G(θ 1 , x) = x -S b 1 a T 1 x (see
r j (1) = (Id IR n -S b 1 P a 1 ) (x j ) r j (2) = (Id IR n -S b 2 P a 2 ) (r j (1)) ... ... r j (d) = (Id IR n -S b d P a d ) (r j (d -1)).
that can be rewritten as:

r j (d) = (Id IR n -S b d P a d )...(Id IR n -S b 2 P a 2 )(Id IR n -S b 1 P a 1 )(x j ) = 1 k=d (Id IR n -S b k P a k ) (x j ) .
Finally, it appears that the representation of X is defined by G(θ, x) = 0 with:

G(θ, x) = 1 k=d (Id IR n -S b k P a k ) (x) , (14)
where θ = (θ 1 , ..., θ d ), the error being r j (d) = G(θ, x j ).

Properties

The model [START_REF] Gifi | Nonlinear multivariate analysis[END_REF] has the following properties : (1) with d axis, G(θ, x) = 0 defines a d-dimensional manifold, (2) the errors are decreasing, (3) with d = n, the model is exact. Proof is given in [START_REF] Girard | Design and statistical learning for nonlinear auto-associative models[END_REF][START_REF] Girard | Position of principal components analysis among auto-associative composite models[END_REF]. Let us note that linear PCA shares these properties. Like for PCA, we define the information ratio associated to the d-dimensional model:

K d = 1 - j r j (d) 2 / j x j 2 .
We deduce from the previous proposition that the K d series are increasing and that K n = 1. K d allows to choose the model dimension for a given information ratio.

Let us come back to the condition [START_REF] Friedman | Projection Pursuit Regression[END_REF]. It justifies the residual errors definition. This definition comes directly from the model equation G(θ, x) = 0 by setting r j (d) = G(θ, x j ). This way of defining an approximation error may be ill-conditioned. For instance when G(θ, x) = 0 defines a quadratic curve, r j does not come from an orthogonal projection on the curve (except in the circle case), and r j can even be infinity for some points close to the curve (see [START_REF] Bookstein | The measurement of biological shape and shape change[END_REF]). However, in our case, for a given axis, we avoid this problem thanks to the orthogonality [START_REF] Girard | Auto-Associative Models and Generalized Principal Component Analysis[END_REF].

Experimental results

First experiment

This is an academic example. Fig. 1(a) shows a sample of a set of curves obtained by translating a given curve. Let us denote this set as {g j (s), j = 1, ..., N } where N = 100 and j corresponds to the translation index. All these curves are sampled on a same interval. For each of them, this leads to a vector x j = ((g j (s 1 ), ..., g j (s n )), with n = 50. This set corresponds to a cloud X in IR n located on a 1-dimensional manifold. The number of examples N may look small compared to the space dimension n, but the important thing here is to have a large number of samples compared to the intrinsic dimension of X . We can imagine what this manifold looks like by projecting it on the first principal plane computed by PCA (Fig. 1(d)). Using the 1dimensional PCA approximation, we get a straight line. This representation is not sufficient. To get a good approximation of the one-dimensional set X , five axes are necessary. This leads to a five-dimensional linear subspace representing 95% of the cloud variance. PCA simulations provide points in IR n are not representative of the training set (Fig. 1(b)), illustrating in this case, the very poor PCA generalization ability. Now, let us consider the NLPCA modeling. We do not know if there is a parametrization axis. The translation, which was used to build the set of examples cannot be written as a linear combination of the coordinates. In order to search for a parametrization axis, the algorithm A ′ was used for the index maximization. The resulting axis a 1 yields the neighborhood preservation of 93 points among the 100 initial points. The generalization criterium ((15) in Appendix) has given an optimal spline smoothing for ν = 29, (Note that the standard cross-validation criterium [START_REF] Girard | Designing nonlinear models for flexible curves[END_REF] gives ν = 18). This generalization criterium was computed using M = 8000 simulated values {ȗ k , k = 1, ..., M } on the axis a 1 , with respect to the probability density of the projection of X on this axis. It gives M points S b (ȗ k ) in IR n located on the manifold approximating X . Fig. 1(d) shows the projection of theses points on the PCA principal plane. The corresponding simulated curves are very close to the ones of X , (Fig. 1(d)). Let us note that a bad choice of the number of knots (ν = 96) leads to a very poor generalization behavior (Fig. 4).

Second experiment

The object of interest is a lamp for which N = 45 appearances have been obtained by varying the azimuth and elevation of the viewpoint1 (Fig. 5). This experiment is close to that of Murase and Nayar [START_REF] Murase | Visual learning and recognition of 3D objects from appearance[END_REF] who have shown the nonlinear nature of such data. A one dimensional NLPCA model has been selected. The algorithm A ′ used with T 0 = 1 and λ = 0.995, has converged after Figure 4: Simulation with a wrong number of knots. 1000 iterations. It tooks 9 mn on a Pentium 233 M hz. The axis a 1 keeps the neighborhood of 33 points among the 45 points of X . In Fig. 6, the representation of the lamp by NLPCA is more accurate than the PCA representation which shows blurred contours and detail removing. In particular the knob of the lamp is not depicted.

Third experiment

X is composed of N = 400 face images from 40 individuals2 (Fig. 7, line 1). Here, the quality of the representation is measured from the approximation error. To represent faces with 20% of mean error, the PCA model requires d = 210 whereas the NLPCA model requires d = 89 for ν = 2 and d = 65 for ν = 5. In this experiment, let us note that PCA and NLPCA axes are close. Fig. 7 (lines 2 and 3) show representation by PCA and NLPCA for d = 89 in the two cases. The mean error is 35% for the PCA representation and 20% for the NLPCA one. As for the lamp experiment, we see that the NLPCA model keeps many details in contrast to the PCA model which yield severely blurred faces.

Conclusion

We proposed a parametric model to approximate a set of points X which are non-linearly distributed in a large multidimensional space.

From the geometric point of view, this representation is done using a d-dimensional manifold, computed by a Projection Pursuit algorithm. The main points of this method are the following : The axes are obtained by optimizing an index which preserves the neighborhood structure of X . When the optimized index is large, the smoothing makes sense because the axis provides a natural order to the points. The second point deals with the ability to control the manifold dimension using the information ratio. Then, the last point is the determination of the smoothing rate using a generalization criterion. These operations are fully automatic, they do not ask the user to adapt some parameters. Finally, when the index and the information values are close to their maximum value, then the model can be considered as valid. smoothing behaviour between the points u j by dropping successively one point u j at a time, as follows :

C(ν) = 1 N N j=1 x j -S b b(ν,j) (u j ) 2 , (16)
where b(ν, j) is the least-square solution based on the data X \ {x j }. However, when X suffers a curse of dimensionality [START_REF] Huber | Projection Pursuit[END_REF], this criterium seems to be less adapted than the generalization error [START_REF] Girard | Building and training radiographic models for flexible object identification from incomplete data[END_REF].

Figure 1 :

 1 Figure 1: (a) Translated curves. (b) Simulation using PCA. (c) Simulation using NLPCA. (d) Projection on the first PCA plane.

Figure 2 :

 2 Figure 2: (a) Gaussian cloud. (b-c-d) Non Gaussian cloud.

( 1 )Figure 3 :

 13 Figure 3: Random walk on the regions.

  choice of S b is much more simpler than the index choice. We express S b as a multivariate spline function S b : IR → IR n [10]. This representation is well-known. The i-th coordinate of S b is noted S (i) b (i) with b = (b (i) ) 1≤i≤n , We approximate the ith coordinate {x(i) 1 , ..., x

  . Go back to 2, if the residual errors are too large, with k ← k + 1; else d = k and end.

Figure 5 :

 5 Figure 5: Several appearances of a lamp.

Figure 6 :

 6 Figure 6: (a) An appearance. (b) PCA representation. (c) NLPCA representation.

Figure 7 :

 7 Figure 7: First line: 3 appearances of a face. Second line: PCA representation. Third line: NLPCA representation.

Centre for Intelligent Systems, University of Plymouth http://www.cis.plym.ac.uk/cis/3Darchive.html.

Olivetti and Oracle Research Laboratory http://www.cam-orl.co.uk/face database.html.

Let us note that the modelling process leads to the estimation of a large number of parameters : axis parameters and spline parameters. However, there is no over-parametrization but on the contrary a dimensionality reduction. This appears clearly when the method is used in a compression framework. Once the model estimated on X , the spline parameters are definitively saved. A new observation compression x only requires the transmission of d scalars, the projections of x on the d axes. The restoration is given by the points on the manifold relative to the d scalars thanks to the model generalization properties. The compression rate is very important.

In fact, our research has been motivated by the need to synthesis artificial observations x. In this framework, the set X is seen as an observation set of a flexible "object" and we are more interested in the simulation of deformation than in the approximation. Precisely, simulating "object" deformations provides useful prior information to solve some inverse problems [START_REF] Girard | Building and training radiographic models for flexible object identification from incomplete data[END_REF].

Appendix

Choice of the number of knots. The problem is a classical one, numerous works have dealt with it. A well-known solution is based on the cross-validation error [START_REF] Whaba | How to smooth curves and surfaces with splines and cross-validation[END_REF]. As it was said above, our final goal is to simulate X in order to get simulations that "look like" the observations X . With the approximation error alone, the model quality is measured only at the u j points. The use of the model can lead to unrealistic simulations away from these points, that is to say, very different from the observations X (see Fig. 2(c)). So, we need a criterion which takes into account the smoothing behavior between the points u j . We propose to use the following generalization error.

Definition 1 Let U be the random variable X projected on the axis a : U = a T X and f (u) its probability density on IR. The theoretical generalization error is

where Q 2 is the square of a distance, E is the mathematical expectation with respect to f (u) and b is the least-square solution [START_REF] Friedman | Exploratory Projection Pursuit[END_REF] .

In practice, an empirical version G emp of G is proposed [START_REF] Vapnik | Estimation of dependences based on empirical data[END_REF]. Let f (t) be an estimation of the density f (t), (see [START_REF] Silverman | Density estimation[END_REF]). In order to define this empirical version, we simulate U according to the f density. Let ȗ be such a simulation and φ(ȗ) the point u j which is the closest neighbor of ȗ on the axis a : φ( t) = arg min j (ȗ -u j ) 2 . The empirical generalization error is then defined using M occurrences (ȗ 1 , ..., ȗM ) of U as follows :

G emp is a quadratic distance between the simulations {S b b(ν) (ȗ k )} of X and the initial data X .

(At this point, we can see how X is simulated using the random variable U ). We choose the number of knots ν that minimizes this criterion. (Let us note that when ν is fixed, the knots are chosen so as to give to each interval defined by two neighbor knots the same probability for f (t)). We could easily show that this criterion can be splitted into a variance term and a bias term and that the ν achieves a trade-off between these two terms [START_REF] Geman | Neural networks and the bias/variance dilemma[END_REF]: when ν is too large, the bias is small (good data approximation) but the variance is large (bad generalization). Finally, let us note that the classical cross-validation criterium C takes also into account the