Hiba Alawieh

Nicolas Wicker

Christophe Biernacki

Projection under pairwise distance control

Keywords: Data visualization, dimensionality reduction, multidimensional scaling, principal component analysis, kernel principal component analysis

Visualization of high dimensional and possibly complex data onto a low-dimensional space is often difficult. Several projection methods have been already proposed to display such high-dimensional structures on a lower-dimensional space, but the information lost is not always considered. Here, a new projection paradigm is presented to describe a non-linear projection method that takes into account the projection quality of each projected point in the reduced space, this quality being directly available at the scale of this reduced space. More specifically, this novel method allows for a straightforward visualization of data in R 2 with a simple reading of the approximation quality, and thus provides a novel variant of dimensionality reduction.

Introduction

Several domains in science use data with large numbers of variables in their studies such as in biology [START_REF] Cheung | Classification approaches for microarray gene expression data analysis[END_REF][START_REF] Golub | Molecular classification of cancer: class discovery and class prediction by gene expression monitoring[END_REF], chemistry [START_REF] Svante | Multivariate Data Analysis in Chemistry[END_REF], geography [START_REF] Van Der Hilst | Seismo-stratigraphy and thermal structure of earth's core-mantle boundary region[END_REF]) and finance [START_REF] Jagannathan | Risk reduction in large portfolios: why imposing the wrong constraints helps[END_REF]. These data can be viewed as a large matrix and extracting results from this type of matrix Hiba Alawieh. Email: alawieh.hiba@gmail.com is often difficult and complicated. In such cases, it is desirable to reduce the number of dimensions of data by conserving as much information as possible from the given initial matrix.

Different types of multivariate data analysis methods have been developed to study these data such as dimensionality reduction, variable selection, cluster analysis and other methods. Typically, dimensionality reduction is used to summarize the data with variable selection used to choose the pertinent variables from the set of candidate variables and cluster analysis used to group the objects or variables. In our study, we focus on dimensionality reduction. Dimensionality reduction techniques can be used in different ways, to solely lower the dimensionality to prepare data for other treatments or for data visualization to provide a simple interpretation of the data in R 2 or R 3 .

Due to the difficulties faced by high dimensional data, many methods for data dimensionality reduction and data visualization have been proposed [START_REF] Chan | A survey on multivariate data visualization in Science and technology[END_REF][START_REF] Chinchilli | Multivariate Data Analysis: Its Methods[END_REF][START_REF] Dempster | An overview of multivariate data analysis[END_REF][START_REF] Keim | Visualization Techniques for Mining Large Databases: A Comparison[END_REF][START_REF] Mardia | Multivariate analysis[END_REF]. Some of the most common methods include principal component analysis (PCA) [START_REF] Jackson | A Users Guide to Principal Components[END_REF], multidimensional scaling (MDS) [START_REF] Togerson | Theory and methods of scaling[END_REF], scatter plot matrix [START_REF] Cleveland | Dynamic Graphics for Statistics[END_REF], parallel coordinates [START_REF] Inselberg | The Plane with Parallel Coordinates[END_REF]) and Sammon's mapping [START_REF] Sammon | A nonlinear mapping for data structure analysis[END_REF]. Scatter plot matrix and parallel coordinates methods are widely used to visualize multidimensional data sets. An issue with PCA and MDS is that as the number of dimensions grows, important multi-dimensional relationships might not be visualized. Moreover, the quality of projection usually assessed by the percentage of variance (PCA case) that is conserved or by the stress factor (MDS case) is a global projection quality measure and does not give information about local quality.

In some projection methods such as PCA, a local measure is defined to indicate the projection quality of each projected point taken individually. This local measure is evaluated by the squared cosine of the angle between the principal space and the vector of the point [START_REF] Jollife | Principal Component Analysis[END_REF]. A good representation in the projected space is hinted by high squared cosine values. This measure is useful in cases of linear projection, which happens in PCA, but cannot be applied in the case of nonlinear projection. Moreover, linear dimensionality reduction misses important nonlinear structure in the data which does not allow to give powerful results in case of nonlinear configurations. Therefore, many methods have been developed to perform nonlinear projections by nonlinearizing a linear dimensionality reduction or by using manifold learning methods.

The nonlinearization of linear dimensionality reduction is applied to extract nonlinear principal components. Kernel PCA is one of the most popular methods in this domain, which integrates a kernel function to determine principal components in different highdimensional space [START_REF] Schölkopf | Nonlinear Component Analysis as a Kernel Eigenvalue Problem[END_REF]. Manifold learning methods are an approach to construct a matrix using the neighborhood information and take a spectral decomposition to find a nonlinear embedding (like Locally Linear Embedding LLE, Isomap algorithm etc) [START_REF] Lee | Nonlinear Dimensionality Reduction[END_REF][START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF][START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF].

In this paper, we propose a new nonlinear projection method that projects the points into a reduced space by using the pairwise distance between pairs of points and by taking into account the projection quality of each point taken individually. Nonlinear projection methods cited in the previous paragraph project the points in a feature space which makes the distances between the projected points hard to be interpreted.

In our method, the distances between projected points are related to the initial distances between points, offering a way to easily interpret the distances observed in the projection plane. This projection leads to a representation of the points as circles with a different radius associated to each point. Henceforth, this method will be referred to as "Projection under pairwise distance control". Furthermore, visualization of data in a reduced space is not the only objective of this method. It can serve as a dimensionality reduction method to reduce the number of variables by minimizing the sum of the radii and to then determine the number of variables that can be kept.

The main contribution of this study is to provide a simple data visualization in R 2 with a straightforward interpretation and to provide a new variant of dimensionality reduction. Firstly, the new projection method is presented in Section 2. In Section 3, the algorithms used in solving the optimization problems related to this method are then illustrated. In Section 4 the application of this method to various real data sets is shown. Finally, the conlusions are drawn in Section 5.

Projection under pairwise distance control

Let us consider n points given by their pairwise distances denoted by d ij for i, j ∈ {1, . . . , n}. The objective is to project these points using distances into a reduced space R q by introducing additional variables, called hereafter radii, that indicate the extent to which the projection of each point is accurate. The local quality is then given by the values of the radii. A good projection quality of point i is indicated by a small radius value denoted by r i . It is important to note that both units of d ij 's and r i 's are identical, thus allowing for a direct comparison.

Before presenting our method, an overview of principal component analysis, Kernel

PCA and multidimensional scaling is given to highlight the significance of our method.

2.1. Overview of certain existing methods: PCA, KPCA and MDS

Principal Component Analysis (PCA)

The PCA method is the most used linear projection technique for data visualization and dimensionality reduction. PCA can be stated as an optimization problem involving the squared Euclidean distances [START_REF] Mardia | Multivariate analysis[END_REF]. This optimization problem is the following:

P PCA :              min A∈Mp×q 1≤i<j≤n |d 2 ij -Ay i -Ay j 2 | s.t. rank(A) = m AA T = I p ,
where y i ∈ R p is the original coordinates vector of point i, d 2 ij is the squared distance for couple (i, j) given by y i -y j 2 and A is the projection matrix of dimension p × q with q being the reduced space dimension. By its nature, PCA cannot take into account nonlinear structures, as it describes the data in terms of a linear subspace. To deal with nonlinearity, Kernel PCA, the reproducing kernel Hilbert space variant of PCA, can be used.

Kernel PCA (KPCA)

The idea behind KPCA is to perform PCA in a feature space denoted by F, obtained by a nonlinear mapping of data from its original space into the feature space F, where the low-dimensional latent structure is hopefully easier to discover [START_REF] Schölkopf | Nonlinear Component Analysis as a Kernel Eigenvalue Problem[END_REF].

The mapping function noted Φ is considered as:

Φ : R p → F Y → Φ(Y ) •
The original data y i is represented in the feature space as a function Φ(y i ) = k(y i , .),

where k(., .) is a positive kernel. Similar to PCA, KPCA is based on finding the first q eigenvectors corresponding to the q largest eigenvalues λ i of the Gram matrix K = (k ij ) ij∈1,...,n , where k ij = k(y i , y j ) = Φ(y i ), Φ(y j ) is a chosen positive kernel. Letting V v , for v = 1, . . . , q, the eigenvectors in the feature space and P Φ(yi) the projection of Φ(y i ) onto the subspace V 1 , . . . , V q . The KPCA problem can be represented as a minimization problem with the following error:

E KPCA : Φ(y) -P Φ(y) 2 2 ,
where

P Φ(y) = q v=1 Φ(y), V v V v •
Furthermore, the most well-known and used measure applied to evaluate the projection quality of points for PCA and KPCA is the squared cosine value. Squared cosine values cannot be interpreted at the same time as the distances in the projection because the cosine values do not have a specific unit. More precisely, the visualization of the projection in the reduced space using PCA and KPCA cannot simply be interpreted in terms of original distances between the points. Indeed, in PCA, the cosine values do not provide a quantitative assessment of the error made when considering the distances between the projected points, even less in KPCA where the projected points are in the feature space so the term "distances" is not related to the distances between the points in the original space.

Multidimensional Scaling (MDS)

As with PCA, Multidimensional scaling (MDS) consists of finding a new data configuration in a reduced space. The main difference between these two methods is that the input data in MDS is in the form of a similarity or dissimilarity matrix, called "proximity", representing the proximity between pairs of objects. MDS are developed where the proximities behave like distances or not respectively [START_REF] Borg | Modern Multidimensional Scaling: Theory and Applications[END_REF]Groenen 2005, Shepard 1962). The key idea of MDS is to perform dimensionality reduction in a way to approximate high-dimensional distances denoted by δ ij the low-dimensional distances d ij , where d ij is equal to the distance between x i and x j , the coordinates of i and j in the reduced space. In his original paper on MDS [START_REF] Kruskal | Nonmetric multidimensional scaling: a numerical method[END_REF]), Kruskal proposed the least-squares loss function denoted by "Stress" as follows

Stress = 1≤i<j≤n (d ij -x i -x j ) 2 1≤i<j≤n d 2 ij •
By minimizing the Stress function, we find the best configuration of (x 1 , . . . , x n ) ∈ R q such that the distances fit to the initial distances.

If we consider n variables as r 1 , . . . , r n ∈ R + , the sum of which bounds the stress function, the optimization problem P MDS can be equivalently rewritten as:

P MDS :                      min x1,...,xn∈R q ,r1,...,rn∈R + n i=1 r i s.t. n i=1 r i ≥ 1 n -1 1≤i<j≤n (d ij -x i -x j ) 2 1≤i<j≤n d 2 ij •
Note that the optimal solution of the MDS problem may not be unique [START_REF] Kruskal | Multidimensional Scaling. Series: Quantitative applications in social Sciences[END_REF].

A criterion to determine the local projection quality has been proposed by Born and Groenen called Stress-per-point (SP P ) [START_REF] Borg | Modern Multidimensional Scaling: Theory and Applications[END_REF]. The SP P of point i is given by:

SP P i = n j=1,j =i (d ij -x i -x j ) 2 n j=1,j =i d 2 ij Stress , with Stress = n 1≤i<j≤n (d ij -x i -x j ) 2 n 1≤i<j≤n d 2 ij •
Again, this is difficult to interpret directly on the projection as a distance error because the projected points are not in the same metric as the initial data.

However, we can observe that the constraint on n i=1 r i can be modified to have a stronger control on each d ij in the following way: |d ij -x i -x j | ≤ r i + r j where x i and x j are the projected coordinates of points i and j.

Therefore, our objective is to propose a new nonlinear projection method that individually controls the projection of points and provides a graphical representation in the same metric as the original space with an error associated to each point.

Our proposal: Projection under pairwise distance control method

Let x 1 , . . . , x n be the coordinates of the projected points in R p and x i -x j the distance between two projected points (i, j). Radii are introduced in this paper to assess how far x i -x j is from the given distance d ij . Indeed, for the couple (i, j), we are aiming for a x i -x j value close to d ij , which should imply a small radius (r i , r j ).

Figure 1 depicts this idea: for each point i ∈ {1, . . . , n}, the projection of i belongs to a sphere with center x i and radius r i such that for each couple (i, j) ∈ {1, . . . , n} we

have x i -x j -(r i + r j ) ≤ d ij ≤ x i -x j + r i + r j .
Radii for uncertainty metric: The idea presented above can be expressed by finding the value of radii that satisfy these two constraints:

• n i=1 r i is minimal. • d ij ∈ [ x i -x j -r i -r j ; x i -x j + r i + r j ], for 1 ≤ i < j ≤ n• x i x j r i r j d ij • • Figure 1.. Example of radii for bounding of the original distance d ij
The projection under pairwise distance control problem can be written as the following optimization problem:

P r,x :        min r1,...,rn∈R + ,x1,...,xn∈R q n i=1 r i s.t |d ij -x i -x j | ≤ r i + r j , for 1 ≤ i < j ≤ n
Linear optimization program using fixed coordinates (x 1 , x 2 , . . . , x n ): Of course, by fixing the coordinates vectors x i for all i ∈ {1, . . . , n} using principal component analysis or any other projection method, the optimization problem can easily be solved in (r 1 , . . . , r n ) using linear programming. This problem can be written as follows:

P r :        min r1,...,rn∈R + n i=1 r i s.t |d ij -x i -x j | ≤ r i + r j , for 1 ≤ i < j ≤ n
It should be noted that a solution for problem P r always exists. Indeed, to satisfy the constraints it is sufficient to increase all r i . Thus, for any method producing points in a reduced space as PCA for instance, we can compute the radii as a post-processing to assess the local quality of the projected points.

P r,x is a non-convex optimization problem: For any dimension p, even with p = 1, note that the optimization problem P r,x is not convex. Indeed, to easily illustrate this fact, we take the function g(x, y) = |d -x -y | considering two solutions (x 1 , y 1 ) = (0, 2) and (x 2 , y 2 ) = (3, 1) with d equal to 2. Thus, we have g(x 1 , y 1 ) = 0 and hyper-space [START_REF] Boggs | Sequential quadratic programming[END_REF]; the active-set method, which is composed of two phases, wherein for the first phase (the feasibility phase) the objective function is ignored while a feasible point is found for the constraints, and in the second phase (the optimality phase) the objective function is minimized while feasibility is maintained [START_REF] Wong | Active-Set Methods for Quadratic Programming[END_REF][START_REF] Cristofari | A Two-Stage Active-Set Algorithm for Bound-Constrained Optimization[END_REF]. The choice of optimization method to use to achieve optimality of the optimization problem is essential and depends on many factors such as the type of problem, desired quality of solution, time limit and availability of the algorithm implementation etc. In fact, all of the methods cited above can be used in optimizing problem P r,x which is a constrained optimization problem having inequality constraints and they are all available in MATLAB using the function "fmicon" for constrained nonlinear optimization problems. Having small radii is the main constraint in our optimization problem, thus the objective is to obtain good solution within a reasonable and practical timeframe. Therefore, a method that balances time and quality of the solution is required.

g(x 2 , y 2 ) = 0 but g x 1 + x 2 2 , y 2 + y 2 2 = d - x 1 + x 2 2 - y 1 + y 2 2 = |2 -0| = 2 which is larger than g(x 1 , x 2 ) + g(y 1 , y 2 ) 2 =
Another strategy of use: Dimensionality reduction One of the main objectives of high-dimensional data studies is to choose, from a large number of variables, those that are important for understanding the underlying studied phenomena. In addition to visualization, our aim can thus be to reduce the dimension rather than to visualize data in R 2 . Therefore, the proposed method can serve to reduce the number of variables by taking into account the value of n i=1 r i . Indeed, by solving the problem P r,x using different dimension values, we can choose the dimension with respect to the local projection quality promoted in this study.

A toy example for illustrating our method

Let us apply the proposed projection method to a simple example by taking a tetrahedron with all pairwise distances equal to 1. For problem P r , the coordinates of points

x i for i = 1, . . . , 4 are obtained using multidimensional scaling. The optimization was carried out using the MATLAB software with the optimization toolbox for linear and nonlinear optimization problem used for problems P r and P r,x , respectively. The value of 4 i=1 r i is equal to 0.7935 for problem P r and 0.4226 for P r,x . It is clear that problem P r,x gives better solutions than problem P r with smaller radii, which indicates better projection quality of points. This result is shown in Figure 2, which depicts the solution obtained using P r and P r,x .

In Figures 2a and2b, the circles with different radii indicate the quality of projection for each point. The circle color is related to the radius value, the shades of gray lie between white and black, the smaller the radius, the darker the circle. The points that have circles with small radii are also considered as projected points. Note that the points represented as points and not as circles are very well projected, having radii almost equal to zero.

In Figure 2b, just one circle appears indicating that the projection quality using problem P r,x is better than when using problem P r . In Figure 2a, half of the points are well projected whereas the other half have large radii, indicating that they are not well projected. Moreover, it is worth noting that the three outer points all have radii equal to 0, which indicates that they are all perfectly placed with respect to one another.

In Figure 2b, the distances between the three points that are very well projected are equal to the distances between these points in their original space (d kl = x k -

x l
where k and l are two very well projected points) whereas the distances from the badly projected points to the perfectly projected points are not yet conserved.

Therefore, using the proposed method, we have succeeded in conserving half of the original distances in the new projection plane and the other half have been changed

to fit the new configuration. If we now apply the proposed method to the distances obtained by MDS to find the radius of each projected point (Figure 2a), it can be noted that one distance is conserved as the original distance and the other five distances are changed which indicates that the proposed method projects the points well by conserving the distances between the points as much as possible.

It is also important to note that, in general, our method is not only a nonlinear projection method with local quality measure, but it can act as a new tool to give the local quality of projection for the classical projection methods using the radii by solving problem P r . It can be used outside our method as post-processing of classical methods.

Connection with existing methods

Multidimensional fitting (MDF) [START_REF] Berge | Multidimensional fitting for multivariate data analysis[END_REF]) is a method that modifies the coordinates of a set of points in order to make the distances calculated on the modified coordinates similar to a given set of distances on the same set of points. The so-called "target matrix", the matrix that contains the point coordinates and "reference matrix"

is the matrix that contains the given distances.

Let us take X = {x 1 | • • • |x n }, the target matrix of coordinates and D = {d ij }, the reference matrix of distances. The objective function of the MDF problem is given by:

1≤i<j≤n |d ij -x i -x j |.
Proposition 2.1. Problem P r,x is bounded from below by

1 n -1 1≤i<j≤n |d ij -x i -x j |
where x 1 , . . . , x n is the optimum for the associated MDF problem.

Proof. By summing all the constraints of problem P r,x , we obtain:

1≤i<j≤n |d ij -x i -x j | ≤ 1≤i<j≤n (r i + r j ) = (n -1) n i=1 r i • So, n i=1 r i ≥ 1 n -1 1≤<i<j≤n |d ij -x i -x j |,
which concludes the proof.

Optimization tools for performing the proposed method

Problem P r,x can be solved using different initialization points for the coordinate matrix X. In this section, we first discuss the different initialization points of the proposed optimization problem and then propose two algorithms to be used in our optimization.

Initialization point for problem P r,x

Different solutions of problem P r,x can be obtained using different initial values of matrix X. We have considered three possibilities:

1-Initial point using a known projection method The first possibility is to use the matrix obtained by PCA or another projection method. The choice of method must be based on the type of data. In this application, we use PCA for quantitative data and MDS for categorical and functional data.

2-Initial point using squared distances

The optimization problem P r,x can be changed by taking the squared distances between points instead of the distances.

Rewriting r 2

i as R i , the problem is changed into

P R,x :        min R1,...,Rn∈R + ,x1,...,xn∈R k n i=1 R i s.t. |d 2 ij -x i -x j 2 | ≤ R i + R j , for 1 ≤ i < j ≤ n.
This transformation is interesting because if the constraints of problem P R,x are satisfied, the constraints of problem P r,x will also be satisfied. Indeed,

|d 2 ij -x i -x j 2 | ≤ R i + R j = r 2 i + r 2 j •
If without loss of generality, d ij ≥ x i -x j , we obtain:

(d ij -x i -x j ) (d ij + x i -x j ) ≤ r 2 i + r 2 j ≤ (r i + r j ) 2 ⇒ |d ij -x i -x j | 2 ≤ (r i + r j ) 2 ⇒ |d ij -x i -x j | ≤ (r i + r j )•
In this way problem P R,x can serve as an initial step in solving problem P r,x .

3-Initial point using an improved solution of problem P r This strategy is more involved. First, we need two properties that provide a way to improve the optimization results of problem P r,x .

Proposition 3.1. Let us consider a point x i such that for an index j, the following inequality is saturated:

|d ij -x i -x j | ≤ r i + r j ,
and the other inequalities involving i are not saturated. The corresponding solution can then be improved by moving x i along the line x j -x i in order to decrease r i and

|d ij -x i -x j |.
Another manner to improve the resolution of problem P r,x is to perform a scale change by multiplying the coordinates x i , for i = 1, . . . , n, by a constant a ∈ R. Thus, the new optimization problem is given by: P r,a :

         min r1,...,rn,a∈R + n i=1 r i s.t. |d ij -a x i -x j | ≤ r i + r j • Proposition 3.2. Let
r i with η(a) = 1≤i<j≤n |d ij -a x i -x j |, then ∃ r1 , . . . , rn a solution of P r,a such that n i=1 ri < n i=1 r i .
The new initial point called X imp , is the improved solution given by using these two properties as follows:

• Firstly, improving the solution of problem P r by solving problem P r,a and using proposition 3.2.

• Secondly, improving the solution of problem P r,a using proposition 3.1.

A deterministic strategy: Algorithm 1

As discussed, three possibilities of coordinate matrix X can be used as the initial point:

1-Coordinates given by PCA or MDS: X PPCA/MDS is the coordinate matrix obtained by applying PCA or MDS and r Pr is a vector that contains the radius of each point obtained by solving P r .

2-Coordinates given by squared distances: X P R,x is the coordinate matrix obtained by solving problem P R,x and R P R,x = r 2 P R,x is a vector that contains the squared radius for each point obtained by solving the subsequent P R,x problem.

3-Coordinates given by improving the solution of problem P r : X imp is the coordinate matrix obtained by improving the previous solution using Proposition 3.1 and r imp is a vector that contains the radius of each point obtained after each iteration of solving problem P r,a

Finding these matrices requires solving the following optimization problems: P r , P R,x and P r,a . Problems P r and P r,a are both constrained linear optimization problems that can be solved using interior-point or simplex algorithms, which are the most widely used algorithms for linear programming. The interior-point algorithm uses a primal-dual predictor-corrector algorithm and the simplex algorithm uses a systematic procedure for generating and testing candidate vertex solutions to a linear program [START_REF] Murty | Linear programming[END_REF]). On the contrary, problem P R,x is a nonlinear optimization problem that can be solved using one of the nonlinear optimization algorithms cited in Section 2.2. All these algorithms are available in MATLAB using the optimization toolbox and can be used for the corresponding problem.

To find the best solution of problem P r,x , we solve it with the three different initialization matrices described above. For this task, we define Algorithm 1 that gives the best solution using the different coordinate matrices. This algorithm consists of two steps, an initialization step and an optimization step. The initialization step offers three different coordinate matrices to be used in the optimization step as an initial point to quickly reach the best solution. During the optimization step, problem P r,x is solved using one of the nonlinear optimization algorithms mentioned in Section 2.2, starting each time with one matrix of the three initial matrices already found.

Thus, for Algorithm 1, described below, the three different initialization matrices are tried and then the best one is chosen that gives the minimum value of n i=1 r i .

Algorithm 1

Input: D: distance matrix, N : number of iterations.

Initialization step

Project the points using PCA or MDS. Solve P r using a linear optimization method. Obtained solution: (X PPCA/MDS , r Pr ). Solve P R,x using a nonlinear optimization method and starting from the solution of P r obtained at the previous step. Obtained solution: (X P R,x , R P R,x ). X imp ← X PR,X . for t = 1 to N do Solve P r,a starting from X imp using a linear optimization method.

Improve the solution of P r,a . Obtained solution: (X imp ,r imp ). end for Optimization step Optimize P r,x using a nonlinear optimization method and starting from X PPCA/MDS , X P R,x and X imp . Choose the minimal solution obtained by these three different starting points.

A stochastic strategy: Algorithm 2

Problem P r,x is a hard problem, thus it is natural to resort to stochastic optimization methods. In the present case, we resort to the Metropolis-Hastings algorithm (Johansen and Evers 2007) which allows us to build a Markov chain with the desired stationary distribution. The challenging parts are the choice of the proposal distribution and the necessity to solve the problem P r at each iteration. Specifically, the

Metropolis-Hastings algorithm requires:

1-A target distribution:

The target distribution is related to the objective function of problem P r,x and is given by:

π(x) ∝ exp -E(x) T
, where E is a function in R given by:

E(x) = n i=1
r i , where {r 1 , . . . , r n } is the solution of problem P r with fixed x.

The variable T is the temperature parameter, to be fixed according to the value range of E.

2-A proposal distribution:

The choice of the proposal distribution is very important to obtain meaningful results. It should be chosen in such a way that the proposal distribution approaches the target distribution. The proposal distribution q(X → .) is constructed as follows, giving priority to the selection of points involved in saturated constraints:

• For each point i, choose a point j (i) with probability equal to:

P j (i) = λ exp -λ(r i + r j (i) -|d ij (i) -x i -x j (i) |) n k=1,k =i λ exp (-λ(r i + r k -|d ik -x i -x k |))
•

• Choose a constant c ij (i) using Gaussian distribution N k (0, σ).

• Generate a matrix X * by moving each vector x i of matrix X t-1 as follows:

-If d ij (i) -x i -x j (i) > 0 then x * i = x i + |c ij (i) |L ij (i) .
-else

x * i = x i -|c ij (i) |L ij (i) ,
where

L ij (i) = x i -x j (i) x i -x j (i) •
3-A linear optimization problem:

For the matrix X generated at each iteration, we solve the linear optimization problem P r and we choose finally the matrix X and the vector of radii which give the smallest value of n i=1 r i .

Algorithm 1 and Algorithm 2 are both implemented in MATLAB and a code for each algorithm can be provided by the authors upon request.

Numerical applications

The projection method presented has been applied to different types of real data sets and also to a simulated data set to illustrate its practical interest.

Experimental setup

In practice, we have tested the proposed method on different simulated and real data sets by solving the optimization problem P r,x using Algorithm 1 in addition to the proposed Metropolis-Hastings algorithm (Algorithm 2). A distance matrix is required each time. For quantitative data, the Euclidean distance between points y i ∈ R p , for i = 1, . . . , n, is computed by the known formula

d ij = p k=1
(y ik -y jk ) 2 . For categorical data, the distance between two points (i, j) is given through the Eskin similarity measure [START_REF] Boriah | Similarity Measures for Categorical Data: A Comparative Evaluation[END_REF]) computed by the formula

p ij = Q t=1 w t p t ij where p t ij =        1 n 2 t n 2 t + 2 if i t =
d ij = s ii -2s ij + s jj .
In addition, to compute the distances between curves of functional data, we have chosen a measure of proximity similar to that studied by [START_REF] Ieva | Multivariate functional clustering for the analysis of ECG curves morphology[END_REF]. In their paper, the authors develop a proper classification designed to distinguish the grouping structures by using a functional k-means clustering procedure with three sorts of distances. For our work we chose one of these three proximity measures as their results are similar. The proximity measure chosen between two curves F i and F j is the following:

d 0 (F i , F j ) = T (F i (t) -F j (t)) 2 dt.
This measure is calculated using the function metric.lp() of the fda.usc package for the R software (Febrero-Bande and Oviedo de la Fuente 2011).

For problems P r and P r,a , we first applied PCA for quantitative data and MDS for categorical and functional data; a linear programming package, called "linprog" which solves linear programming problems, was then used to solve the optimization problems with an interior-point algorithm. Problems P r,x and P R,x are nonlinear optimization problems; therefore, we used a nonlinear programming package, called "fmincon" which finds minimum of constrained nonlinear multi-variable function, to solve them. The algorithms cited in Section 2.2 can be used here, but we recommend to use the activeset algorithm. Algorithm 2 can provide a good solution if the parameters λ, σ and T are chosen adequately. For instance, λ should be such that the points belonging to unsaturated constraints are chosen with small probabilities. Therefore, we took it equal to 100. For the other parameters σ and T , we took their values in the range from 0.01 to 100. The choice of these numbers is taken after trying different values of σ and T in order to have the best solution that gives a minimal value of

n i=1 r i .
Moreover, the visualization of the projection of each point i in R 2 is represented as a circle having x i as the center and r i as the radius in a two-dimensional space, where the horizontal and vertical axes correspond to the first and the second dimension associated to the projection in R 2 , respectively. The projected point belongs to this circle and this is the specificity of our method. For each data set, the circles obtained for each point after solving the optimization problem P r,x are shown. To compare the projection quality of our representation with that obtained by PCA and KPCA, we used the squared cosine values as projection quality, and for MDS, the Stress-perpoint (SP P ) [START_REF] Borg | Modern Multidimensional Scaling: Theory and Applications[END_REF]. Indeed, for PCA and KPCA, we plotted the projected points indexed by their squared cosine values and for MDS, we used the smacof package in R to compute the stress-per-point and to plot the bubble plot represented the stress-per-point.

A simulation study

To evaluate the performance of projection under pairwise distance control method, we conducted a simulation study. We generated 100 random samples of y i from a 5dimensional multivariate normal distribution with mean 0 and covariance matrix I, the identity matrix, and we calculated the Euclidean distances between pairs (y i ,y j )

for 1 ≤ i < j ≤ n. The projection result was compared with those obtained by KPCA. Figure 3 shows the results of the projection of the simulated data using the proposed method and KPCA. By comparing Figure 3a and Figure 3b, it can be shown that the projection quality of points using KPCA is somehow dependent on the position of the points in the reduced space. Indeed, the projection is likely to give better local projection quality if the projected point is located near to the center (0, 0). On the contrary, this is less visible for the proposed method.

This simulated data illustrates the originality and the efficiency of the proposed method in giving a good local projection quality.

Introducing the real data sets

Four real data sets were used and divided into three categories:

• Quantitative data: Iris and car data sets.

• Categorical data: Soybean data set.

• Functional data: Coffee data set.

The Iris data set [START_REF] Anderson | The Irises of the Gaspé Peninsula[END_REF]) is a famous data set and is presented to show that the projection works as expected. This data set contains 3 classes of 50 instances each,

where each class refers to a species of Irises. The four variables studied in this data set are: sepal length, sepal width, petal length and petal width (in cm). The car data set [START_REF] Saporta | Probabilités, analyse des données et statistique[END_REF]) is a data set studied in the book by Saporta (Table 17.1,page 428). This data set describes 18 cars according to various variables (cylinders, power, length, width, weight and speed).

The soybean data set [START_REF] Stepp | Conjunctive conceptual clustering[END_REF] The coffee data set is a time series data set used in chemometrics to classify food types. It is a functional data set where 56 samples of coffee are available with 286 timestamps for each sample (as a result of spectroscopic analysis). This kind of time series is common in many applications in food safety and quality assurance and was taken from the UCR time Series Classification and Clustering website (Chen et al.

2015)

. Coffea Arabica and Coffea Canephora variant Robusta are the two species of coffee bean that have acquired a worldwide economic importance, and many methods have been developed to discriminate between these two species by chemical analysis [START_REF] Briandet | Discrimination of arabica and robusta in instant coffee by fourier transform infrared spectroscopy and chemometrics[END_REF]. The optimization results for these four data sets are given in Table 1. For each data, the sum of radii n i=1 r i obtained using Algorithm 1 and Algorithm 2 is provided. 1, the solutions of Algorithm 2 for the different data sets are shown to be very close to those obtained using Algorithm 1. Thus, the radii obtained are estimated to be close to the solution of Algorithm 1. Moreover, it is interesting to note here that the number of iterations N in Algorithm 1 has an important role in finding the minimal value of n i=1 r i for problem P r,a and then for problem P r,x and also to reduce the computing time. In fact, the important decrease in the value of n i=1 r i occurred in the first 500 iterations through of 1000 iterations, and then a small decrease occurred after 500 iterations. This small decrease in value of n i=1 r i after 500 iterations shows that a size of 500 iterations can be a good choice for the Algorithm 1 since all the studied data sets are concerned. Indeed, this result can be observed for all data sets presented in our application with approximately 500 iterations.

Iris data set: Figure 4 depicts the result of projection under pairwise distance control for the Iris data set. In the projection of the Iris data set shown in Figure 4, it is interesting to note that two areas are well separated. This corresponds to the wellknown fact that Iris versicolor and virginica are close whereas the species Iris setosa are more distant.

Referring to the original data, the Iris data set contains three classes corresponding Figure 4.. Projection of the Iris data set using projection under pairwise distance control method. Two well separated groups can be observed. The points that have circles with small radii are considered as well projected points.

to the three types of Iris plants and one class is linearly separable from the other two classes. This result clearly appears in our projection.

Moreover, we have compared the local projection quality of PCA, KPCA and MDS with the local projection quality obtained using projection under pairwise distance control. By comparing the projection of PCA with the projection of our method for the Iris data set given respectively in Figures 5 and4, we can say that our method projected the points without giving any importance to any group. Figure 5 depicts a group with small values of the quality measure and another group with high values of quality measure, whereas the radii obtained by projection under pairwise distance control method are distributed in an equivalent way.

For KPCA, we plotted the squared cosine values as circles to make the representation clearer, especially for the Iris data set as the Iris setosa species are projected next to each other. From Figure 6a, we can conclude that in each category, the points that have close quality values are located side by side.

Furthermore, by comparing the proposed projection method with the one obtained by MDS, it can be concluded that, as is the case when using PCA, the points in Figure 6b are projected by giving more importance to the Iris setosa group. Indeed, almost all the red circles (indicating a very good projection) are assigned to the Iris setosa species. Moreover, the comparison of the position of points in the reduced space in terms of distance between points cannot be viewed in this classical method as the points in the reduced space are not in a metric compatible to the initial distances, whereas in our method we have conserved the metric of the initial distances.

Cars data set: The projection of points using projection under pairwise distance control for the car data set is shown in Figure 7. The expensive cars, the "Audi 100", "Alfetta-1.66", "Datsun-200L" and "Renault 30" are well-separated from the lowstandard cars, the "Lada-1300", "Toyota Corolla", "Citroen GS Club" and "Simca 1300". Moreover, we can assert that the projected points obtained using projection under pairwise distance control are well separated as there are no circle intersections.

By comparing our projection with the projection obtained using PCA presented in Figure 8, it can be shown that in the projection of PCA, there is a group with small values of the quality measure located at the center, which corresponds to the cars:

Lanca-Beta, Mazda, Fiat, Simcs and Rancho, and a group with high values of quality measure located far from the center.

Regarding KPCA, we can see in Figure 9a that the points with navy circles are almost all located almost around the same y-axis coordinates and the same applies for the It can also be noticed that the cars Princess, Mazda, Fiat and Peugeot located in the same area with small circles. Therefore, the only conclusion that we arrive at is in relation to the size of the circles and to the quality of the projected points. However, it is not possible to conclude anything about the closeness of these 4 points as the distances here are in the feature space and are not related to the original space.

In Figure 7, we can however conclude that the two cars, the Mazda and Fiat, are well projected in the reduced space, and they have similar characteristics as these two cars are close. The same conclusion can be made for the Peugeot and Princess cars. From this, it is possible to conclude that there is a large difference between the two cars, the "Toyota" and "Renault 3" as the distance between these two cars is significant. Conversely, the distance between the "Lada1300" and "Citroen" is small, thus indicating the closeness of these two cars. Note that these two cars are very well projected, resulting in a very good interpretation of the distance between them.

Therefore, radii are meaningful in our method and give an interpretation about the distances between points whereas the distances between the projected points using PCA, KPCA and MDS are not interpretable. This is a particular strength of our method. Projection under pairwise distance control suggests an absolute interpretation whereas the other methods provide a relative one.

For the qualitative and functional data sets and using MDS, recall the definition of the Gram matrix called B which is equal to X X where X is the coordinate matrix in the reduced space. Thus, it is necessary to verify that the matrix B obtained by the MDS method is semi-definite positive to use the squared cosine as the quality measure because the starting point of optimization is obtained from MDS. After this, in case of positiveness of matrix B, the quality measure can be calculated.

Soybean data set: In the projection of the soybean data set, four classes are shown in Figure 10 and each class contains the disease number of the class. The whole set of points can however be divided in two large classes. Indeed, it is clear that Class 2 is well separated from the other classes as there is no intersection between the circles of Class 2 and the circles of other classes. Moreover, Class 1 can be considered as well separated class from Classes 3 and 4 if the largest circle D 3 is not taken into account.

Classes 3 and 4 are not well separated at all, as there are different intersections between the circles of these two classes. This result is shown in [START_REF] Stepp | Conjunctive conceptual clustering[END_REF] which labels the first two classes as "normal" and the latter two classes as "irrelevant". A comparison of results from projection under pairwise distance control with PCA and KPCA is not possible for this data set because the matrix B is not semi-definite positive. Regarding

Figure 11, it is clear that Class 4 exhibits the worst projection quality, whereas Classes 1 and 2 show better projection quality. Therefore, it is possible to draw the same conclusion for the Iris and car data sets when using MDS as a projection method, the projection quality of points is dependent on the class of the points. In Figure 12, we show that we have succeeded in differentiating the Arabica from Robusta coffee. These two classes are clearly presented, the first class indexed by number 1, corresponding to Arabica coffee, and the second one indexed by number 2, corresponding to Robusta coffee. These classes are not well separated in comparison with the results of quantitative data, since there are many intersections. Therefore, the representation of the points as circles and not as simple points provides more information about the real point classes and shows the points that are at risk of being misplaced in a particular class.

Figures 13a and13b show the projection quality using PCA and MDS respectively.

As all the eigenvalues of matrix B are positive, we can compute the quality measure given by PCA. Comparing the projection quality of PCA and projection under pairwise distance control provided by Figures 13a and12, respectively, it can be seen that the quality of projection of the set of points is quite steady.

Additionally, Algorithm 2 was applied to these data sets. The trace plots of the optimization problem P r,x are shown in Figure 14 after 5000 iterations. It is important to note that the value of the sum of radii n i=1 r i decreases rapidly in the first iterations and stays roughly constant after 1000 iterations for the different data sets, with the Finally, the computer speed time of the proposed method is compared with that using the classical methods. Table 2 shows the computer speed time for the four data sets using PCA, KPCA, MDS, Algorithm 1 and Algorithm 2. It is clear that our method takes more time than the existing methods. However, Algorithms 1 and 2 are expected to significantly increased by using the C++ programming language (instead

Dimensionality reduction results

Our method can also be directly used to reduce the dimensionality of data (possibly using it beyond visualization in R 2 ). This only requires solving problem P r,x using different dimension values. In Figure 15, the values of n i=1 r i were plotted as a guide for choosing the reduced number of variables. This figure shows the values of n i=1 r i for the different data sets using different dimensions. It is clear that the value of n i=1 r i decreases when the dimension increases. Indeed, the sum of radii n i=1 r i decreased rapidly in low dimensions and then decreased slowly when the dimension increases.

The main problem, which is widely posed in dimensionality reduction methods, is the determination of the number of components that need to be kept. Many methods have been discussed in the literature [START_REF] Besse | PCA stability and choice of dimensionality[END_REF][START_REF] Jollife | Principal Component Analysis[END_REF]) to determine the dimension of the reduced space, relying on different strategies related to a good explanation or a good prediction. Thus, with our method the choice of the reduced space dimension is related to the local projection quality of points and how much the user is interested in the projection quality of points. Regarding the quantitative data sets (Iris and car), if the main objective of the user is to obtain a very good projection quality, then a choice of three components against four for Iris data set, and six for the car data set can be a good choice, as the value of n i=1 r i is small and there is not a large difference between this value and the values for higher dimensions. For the coffee data set, a dimensionality reduction from 56 sample time series down to 6 simple extracted features is considered as a good choice. As for the soybean data set, a reduced space dimension equal to 4 dimensions can be considered as an appropriate reduced space.

A comparison of our results with the existing results shows a coherence between them. For the Iris data set, [START_REF] Chiu | Method and Software for Extracting Fuzzy Classification Rules by Subtractive Clustering[END_REF] and [START_REF] Liu | Chi2: feature selection and discretization of numeric attributes[END_REF] concluded that the number of variables could be reduced to 2 as the petal length and petal width variables are the most important variables from all the variables. For the car data set, Saporta (2006) (Table 7.4.1 page 178) noticed that the conservation of two dimensions led to the explanation of 88% of inertia, where the inertia term reflects the importance of a component. Theredore, these results seem very similar to our results, with the important decrease located between dimensions 1 and 2. The other reductions are negligible for these two data sets. A selection of variables was studied on time series coffee data set by [START_REF] Andrews | Variable Selection for Clustering and Classification[END_REF]. Using several analysis methods, the number of selected variables ranged between 2 and 13. This result is also seen using our method, a number of reduced variables taken between 2 and 9 gives a good projection. Regarding the soybean data set, Dela Cruz shows in his paper Dela [START_REF] Cruz | Comparative Study of Data Mining Classification Techniques over Soybean Disease by Implementing PCA-GA[END_REF] that the 35 attributes can be reduced to 15. With our method, we have succeeded in reducing the attributes to 6 by having a very good projection of points.

Hence, the results presented confirm that the dimension nonlinearly can be reduced while assessing a reasonable number of dimensions at the same time.

Conclusion

The purpose of this paper was to outline a new nonlinear projection method based on a new local measure of projection quality. Of course, in some projection methods, a local measure is given but this measure cannot be applied unless in cases of linear projections, and even then it is not efficiency for graphical representation.

The quality of projection is given here by additional variables called radii, which enable bound on the original distances to be obtained. We have also shown that the idea can be written as an optimization problem in order to minimize the sum of the radii under some constraints. As the solution of this problem cannot be obtained exactly, we developed a stochastic optimization method.

This method has several advantages. Firstly, it is a nonlinear projection method that takes into account the projection quality of each point individually. Secondly, the distances between projected points are related to the initial distances between points offering a way to easily interpret the distances observed in the projection plane. The projection quality of each point can even then be used outside our method, as a postprocessing of PCA or MDS for example. Finally, it appears to be efficient in terms of dimensionality reduction for the selection of the dimension of the reduced space based on the local quality of projection.

As perspectives, a lower bound for the optimization problem is needed and this radii approach could also be applied to other methods.

  0 proving non convexity associated to this sample design. Many methods available in the literature propose different ways to solve such optimization problems. Examples include: trust-region-reflective (Conn et al. 2000), which chooses and computes an approximation of the objective function, and then chooses and modifies the trust region and finally solves the trust-region subproblem; sequential quadratic programming (SQP) which solves the optimization problem by addressing a sequence of quadratic programming problems where the Lagrangian function is approximated by a quadratic function and the constraints are approximated by a linear

Figure 2

 2 Figure 2.. Projected points after solving problem P r and problem P r,x . The x-axis and y-axis are dimension 1 and dimension 2, respectively. (a) and (b) show the projection obtained from the solution of problem P r using MDS and of problem P r,x , respectively.

Figure 3

 3 Figure 3.. Projection of the simulated data using the proposed method (a) and Kernel PCA (b). The points that have circles with small radii are considered as well projected points.

  from UCI Machine Learning Repository characterizes 47 soybean disease case histories defined over 35 attributes. Each observation is identified by one of the 4 diseases: Diaporthe Stem Canker (D1), Charcoal Rot (D2), Rhizoctonia Root Rot (D3) and Phytophthora Rot (D4).

4. 4 .

 4 Results from the real data sets 4.4.1. Data visualization in R 2

Figure 5 .

 5 Figure 5.. Projection of the Iris data set using PCA. The values of local projection quality are given for each instance. The values of projection quality for Iris setosa species (Se) vary between 0.97 and 1 indicating then a very good projection quality whereas for Iris versicolor (VeCo) and Iris virginica (ViCa) species, the values of projection quality vary between 0.2 and 0.99 indicating a very large variability in the projection quality.

Figure 6

 6 Figure 6.. Projection of the Iris data set using KPCA (a) and MDS (b). The color convention is as follows: the darker the red color of a particular disk, the better the projection. Inversely, the darker the blue color of a particular disk, the worse the projection.

Figure 7 .

 7 Figure 7.. Projection of the car data set using projection under pairwise distance control.

Figure 8 .

 8 Figure 8.. Projection of the cars data set using PCA. The values of local projection quality are given for each car.

Figure 9

 9 Figure 9.. Projection of the car data set using KPCA (a) and MDS (b).

Figure 10 .

 10 Figure 10.. Projection under pairwise distance control for the soybean data set. Four groups are presented, indexed by D1, D2, D3 and D4.

Figure 11 .

 11 Figure 11.. MDS for the soybean data set. Four groups are presented, indexed by D1, D2, D3 and D4.

Figure 12

 12 Figure 12.. Projection of the coffee data set using projection under pairwise distance control. Two clusters, indexed 1 and 2, indicate the Arabica and Robusta classes respectively.

Figure 13

 13 Figure 13.. Projection of coffee data set using PCA and MDS.

Figure 15

 15 Figure 15.. Scree plots of n i=1 r i for different dimensions for the four data sets.

  r 1 , . . . , r n ; x 1 , . . . , x n be a feasible solution of P r,x , if ∃a such

	n
	that η(a) <
	i=1

Table 1 .

 1 . Solution of problem P r,x for data sets using Algorithm 1 and Algorithm 2.

		n
		r i
		i=1
		Algorithm 1 Algorithm 2
	Iris	16.19	17.2
	Cars	3.27	3.35
	Soybean	3.98	3.93
	Coffee	21.68	21.97
	Based on Table

Appendix

Proof of proposition 3.1

Let us consider a point x i such that for an index j, the following inequality is saturated:

and the other inequalities involving i are not saturated. Then, the corresponding solution can be improved by moving x i along the line x j -x i in order to decrease r i and

Proof. The above condition means that x i is rewritten as

and we look for a such that

Let us now consider the other inequalities corresponding to index pairs (i, k) with k = j. For each of them, ∃a ∈ [a k , a k ] with a k < 0 and a k > 0 such that

as these constraints are unsaturated. Finally, taking a different from 0 in [a , a ] with a = max k a k and a = min k a k , all constraints involving i get unsaturated so that r i can be decreased, thereby decreasing the objective function. Depending on whether a must be negative or positive, we take a = a or a = a respectively.

Proof of proposition 3.2

Let r 1 , . . . , r n ; x 1 , . . . , x n be a feasible solution of P r,x , if ∃a such that η(a) <

Proof. Let us consider r 1 , . . . , r n ; x 1 , . . . , x n a feasible solution of problem P r,x and a, r1 , r2 , . . . , rn a solution of P r,a where a is kept constant. For the solution of P r,a , for each point i, we have a certain saturated constraint associated to point k denoted by C ik(i) , otherwise we can easily saturate it using proposition 3.1. Thus, we have:

Then, |d ik(i) -a x i -x k(i) | = ri + rk(i) ≥ ri . By summing for all points i, for i = 1, . . . , n, we obtain: