Arthur Guillon

Marie-Jeanne Lesot

Christophe Marsala

A Proximal Framework for Fuzzy Subspace Clustering

Keywords: fuzzy clustering, fuzzy subspace clustering, proximal descent

This paper proposes a fuzzy partitioning subspace clustering algorithm that minimizes a variant of the FCM cost function with a weighted Euclidean distance and a non-differentiable penalty term. The form of the cost function suggests to split the optimization problem, taking advantage of the framework of proximal optimization. The expression of the proximal operator for the penalty term is derived and implemented in a new algorithm, PFSCM, which combines proximal descent and alternate optimization. A discussion on the extension of this work to produce sparse estimations of the subspaces is conducted. Experimental results show the relevance of the proposed approach.

Introduction

Cluster analysis is a central task in unsupervised machine learning that aims at partitioning data into groups with strong internal similarity and external dissimilarity. Depending on the considered framework and context, clusters partition the original data set into groups sharing some common property with varying definitions. Most often, cluster analysis informs the user about the prevailing tendencies of the data set.

Subspace clustering [START_REF] Agrawal | Automatic subspace clustering of high dimensional data for data mining applications[END_REF] is a generalization of this task which not only searches for a good partition of the data set, but also identifies the subspaces in which these clusters can be formed, that is, the subspaces in which the points of a cluster are similar, excluding the dimensions or features in which they are dissimilar. The identified subspaces are required to be minimal, but sufficient to describe the clusters they contain.

The identifications of the clusters and their subspaces must be simultaneous: indeed, if either the clusters or their subspaces are known beforehand, the problem reduces to finding the subspaces or correct description of the clusters, or to standard clustering, respectively. In addition, as opposed to feature selection, different clusters are most of the time discovered in different subspaces: the learned metric or similarity measure is local and may differ for each cluster.

As detailed in Section 2 there exist several families of techniques and algorithms to solve the subspace clustering problem, as well as various representations of the subspaces, depending on the intended application of the subspace clustering. This paper belongs to the partitioning paradigm in a fuzzy setting. It produces clusters identified by a center. It focuses on discovering axes-parallel subspaces, which are thus identified by weights on the original data features.

Subspace clustering searches for appropriate description of the clusters using original features, motivating the identification of sparse descriptions of these subspaces, which use as few features as possible. The study of sparse models in subspace clustering is the topic of active research, as detailed in Section 2. 3. In this work, an original cost function formalising these concepts is presented. It adds, to a FCM cost function with weighted Euclidean distance [START_REF] Bezdek | Pattern Recognition with Fuzzy Objective Function Algorithms[END_REF], a penalty term expressing constraints to identify the relevant subspaces. As it is not differentiable, standard optimization techniques such as alternate optimization are not available, and this penalty term requires a new descent procedure which will allow for new constraints on the solutions, such as sparsity. Tools from the proximal descent theory [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF] are thus adapted to the subspace clustering problem. The utilisation of such techniques is still relatively new in machine learning and in clustering in particular [START_REF] Parikh | Proximal algorithms[END_REF].

This paper introduces a new algorithm, called PFSCM, standing for Proximal Fuzzy Subspace C-Means, which uses this framework to solve the subspace clustering problem through the combination of proximal descent and alternate optimization.

This paper is structured as follows: Section 2 sketches related works and the scientific context of subspace clustering. The proposed cost function is presented and studied in Section 3. In Section 4, the implementation of proximal descent is studied to optimize the proposed function, leading to the update equations from which the PFSCM algorithm is derived. Section 5 discusses variants of the proposed cost function based on the introduction of a sparsity-inducing penalty term. PFSCM is then experimentally validated in Section 6. Conclusion and future works are discussed in Section 7.

Related Works

This section first introduces the problem of subspace clustering using standard notations: the data set is denoted {x i , i = 1..n} with ∀i, x i ∈ R d and the cluster set {C r , r = 1..c}. In the following, the related works in fuzzy subspace clustering and sparse methods for subspace clusterings are recalled. These works are presented by order of relevance rather than in chronological order.

Problem Statement

Subspace clustering [START_REF] Agrawal | Automatic subspace clustering of high dimensional data for data mining applications[END_REF][START_REF] Vidal | A tutorial on subspace clustering[END_REF] is an unsupervised task that aims at identifying groups with strong internal similarity while simultaneously learning this similarity, i.e. the set of features or subspace shared by points of a same group. This similarity is local: it represents the subspace in which the cluster points exist. Depending on the field of application, various types of subspaces can be of interest, as hyper-rectangles [START_REF] Agrawal | Automatic subspace clustering of high dimensional data for data mining applications[END_REF], vector subspaces [START_REF] Vidal | A tutorial on subspace clustering[END_REF][START_REF] Aggarwal | Fast algorithms for projected clustering[END_REF] or hyperplanes of low dimension [START_REF] Wang | K-subspace clustering[END_REF] to name a few. Most of the time, subspaces are defined by projections or linear combination of the original axes. This second point of view presents subspace clustering as a generalization of Principal Component Analysis (PCA) in the sense that, in this case, the PCA is local to each cluster.

c 1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • c 2 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • x y z
In this work, subspaces of interest are defined by weights (w rp ) p∈{x,y,z} of the original features according to their importance in the resulting cluster, corresponding to axes-parallel subspaces. Figure 1 illustrates this choice: points of cluster c 1 all exist on the same plane y = 0, which means that the corresponding feature is important to identify the cluster and that it should receive a high weight (w 1y 0). Points of c 1 are also closer along the z axis than they are along the x axis, and thus w 1z > w 1x .

Fuzzy Subspace Clustering

Fuzzy subspace clustering generalizes standard subspace clustering by allowing each data point to be assigned to several clusters with various degrees. It usually relies on the introduction, in standard fuzzy clustering algorithms such as the fuzzy c-means algorithm (FCM) [START_REF] Bezdek | Pattern Recognition with Fuzzy Objective Function Algorithms[END_REF], of subspaces (not necessarily axes-parallel) and weighted distances. Deng et al. [START_REF] Deng | A survey on soft subspace clustering[END_REF] propose a survey on soft subspace clustering models and algorithms.

Keller and Klawonn [START_REF] Keller | Fuzzy clustering with weighting of data variables[END_REF] adapt the FCM cost function by introducing a weight w rp for each cluster C r and dimension p. Denoting (u ri ) ∈ [0, 1] for i ∈ {1, . . . , n} and r ∈ {1, . . . , c} the fuzzy membership degree of x i to cluster C r , c r ∈ R d the center of cluster C r and (w rp ) ∈ [0, 1] the weight of dimension p for cluster C r , they study the following cost function:

J K&K (C, U, W ) = c r=1 n i=1 u m ri d p=1 w v rp (x ip -c rp ) 2 (1)
where m, v ∈ R are fuzzifiers tuned by the user to specify the level of fuzziness of the corresponding parameters and C, U and W are the matrices containing the centers (c r ), the memberships (u ri ) and the weights (w rp ), respectively. The function is minimized under the three constraints:

• (C1) ∀i ∈ {1, . . . , n}, c r=1 u ri = 1 • (C2) ∀r ∈ {1, . . . , c}, n i=1 u ri > 0;
• (C3) ∀r ∈ {1, . . . , c},

d p=1 w rp = 1.
The first two constraints (C1) and (C2) are similar to the FCM ones and ensure that each cluster has an equal importance in the computation of the solution and that it is not trivially empty. Constraint (C3) on the weights (w rp ) is specific to the subspace clustering problem and prevents the trivial solution such that ∀r, ∀p, w rp = 0.

The model used in [START_REF] Keller | Fuzzy clustering with weighting of data variables[END_REF] is a simple adaptation of the FCM cost function to the subspace clustering problem. The Gustafson-Kessel model [START_REF] Gustafson | Fuzzy clustering with a fuzzy covariance matrix[END_REF] is more complex, based on the Mahalanobis distance, capable of identifying non-axesparallel subspaces. Borgelt and Kruse [START_REF] Borgelt | Shape and size regularization in expectation maximization and fuzzy clustering[END_REF] improve upon this model to add shape and size regularization for the clusters.

Sparsity Constraints for Subspace Clustering

The downside of this approach is that each feature typically receives a nonzero weight, even if it plays a small role in the cluster. This not only clutters up the description of the subspaces, but also introduces noise in the cluster analysis of the data, namely the computed centers and memberships. Indeed, as shown by the update equations 3 and 4 in Section 3.2, the weigthing of the dimensions have an influence on the other variables, when they should be identified as irrelevant and not taken into account.

As a good subspace clustering solution consists in minimal yet sufficient descriptions of the clusters and their subspaces, subspace clustering is related to sparsity-inducing techniques. Over the last years, several approaches to subspace clustering heavily relied on sparse solutions to optimization problems in order to identify the correct subspaces. Most of these approaches are related to compressed sensing and are not directly usable in the framework presented in the previous section.

Most notably, Elhamifar and Vidal's Sparse Subspace Clustering (SSC) [START_REF] Elhamifar | Sparse subspace clustering[END_REF] solves an 1 -based optimization problem for each point x i , finding a minimal representation of x i based on its neighbor points. In its simplest setting, this framework makes the assumption that the clusters live in a union of vector subspaces, and thus can be recovered through these local representations (see e.g. [START_REF] Patel | Latent space sparse subspace clustering[END_REF] for a theoretical analysis). A similarity relation is learnt and used to identify the clusters.

This framework has been extended by various authors. Li and Vidal [START_REF] Li | Structured sparse subspace clustering: A unified optimization framework[END_REF] propose a structured approach to SSC, where the sparse representation and the resulting segmentation are refined iteratively through alternate optimization. Instead of looking for sparse representations, Liu et al. [START_REF] Liu | Robust subspace segmentation by low-rank representation[END_REF][START_REF] Liu | Latent low-rank representation for subspace segmentation and feature extraction[END_REF] search for lowrank representations (LRR) of the data points, using the nuclear norm [START_REF] Meyer | Matrix analysis and applied linear algebra[END_REF] as an approximation of the rank function. The underlying idea is that SSC identifies sparse local representations of the data but may have more trouble identifying the global structure of the data in presence of noise, whereas LRR constrains the solution globally.

Closer to the paradigm studied in the present paper, Witten and Tibshirani [START_REF] Witten | A framework for feature selection in clustering[END_REF] reformulate the k-means original problem into a maximization problem with a weighted distance. This allows them to add an 1 -based constraint to the problem in order to produce sparse weight vectors and identify the subspaces. Qiu et al. [START_REF] Qiu | A sparse fuzzy c-means algorithm based on sparse clustering framework[END_REF] extend this framework to fuzzy clustering and compare it to some usual subspace clustering algorithms. Borgelt [START_REF] Borgelt | Fuzzy subspace clustering[END_REF] proposes alternatives to Keller and Klawonn and Gustafson-Kessel algorithms with sparsity inducing fuzzifier functions for u ri and w rp , the second one leading to the selection of principal features or axes. However, both of these contributions cannot introduce sparsity inducing terms without changing the original subspace clustering model.

Proximal Splitting for Fuzzy Subspace Clustering

This section presents a new cost function to model the subspace clustering problem and studies its properties. It is based on the Keller and Klawonn model but moves the constraint (C3) inside the cost function in order to change the optimization dynamic.

A Model for Fuzzy Subspace Clustering

Using the same notations as in Section 2.2, we introduce the following cost function:

J(C, U, W ) = c r=1 n i=1 u m ri d p=1 w 2 rp (x ip -c rp ) 2 F (C,U,W ) +γ c r=1 | d p=1 (w rp ) -1| G(W ) (2)
which is optimized under the classic FCM constraints (C1) and (C2). This function is of the form J(C, U, W ) = F (C, U, W )+γG(W ). As shown in the next section, this form plays an important role in the optimization strategy of J.

Function F is the same as Keller and Klawonn's model, except that the weight fuzzifier v is set to 2 in order to simplify further mathematical analysis. The general case is left to be considered in future works. As it is a special case of Keller and Klawonn's subspace clustering model, function F computes solutions that are close to this model (although the use of a different descent scheme introduces some differences, see Section 6) and thus the function defined in Equation ( 2) aims at preserving these solutions while using a second function G to enforce some additional constraints.

In this paper, function G adds a penalty to the function which prevents the trivial solution W = 0. It can be understood as an inlined version of the constraint (C3), with two differences:

• the hyperparameter γ ∈ R is not a Lagrange multiplier, but is rather set by the user beforehand. Its purpose is to balance out the two terms, and, if large enough, it ensures that the weights W sum up to 1;

• function G is not differentiable and requires new optimization techniques presented in the next section. The advantage is that, as discussed in Section 5, they will allow for more advanced penalties, e.g. for enforcing new properties such as sparsity.

The cost function J in Equation ( 2) thus conveys the idea of finding a solution to the subspace clustering problem (as expressed by F ) with a relaxed constraint (as expressed by G), inspired by 1 -regularization [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF].

Minimization of the Cost Function

As in Section 2.2, an algorithm solving the subspace clustering problem has to find parameters (C * , U * , W * ) minimizing J. As the proposed J is not differentiable, this algorithm cannot be derived from the standard approach of constrained optimization. Instead, the approach we propose is different, and benefits from two observations:

• the second term c r=1 | d p=1 (w rp ) -1| only involves the variable W and thus the optimization of variables C and U is not impacted by the nondifferentiability;

• F is still differentiable in all three parameters and its gradient as a function of W is Lipschitz-continuous: these properties guarantee good performance of well-known optimization algorithms, such as gradient descent.

For the variables C and U , the optimization of F can be done through standard alternate optimization, as for fuzzy c-means: the two update equations for membership degree and cluster centers are derived,

u ri = d 2 1-m ri c s=1 d 2 1-m si where d 2 ri = d p=1 w 2 rp (x ip -c rp ) 2 (3)
and

c rp = n i=1 u m ri • x ip n i=1 u m ri (4)

Proximal Splitting for Weights Optimization

The function J is not differentiable in W on its domain of definition and thus cannot be optimized by standard alternate optimization. In this section, an alternative optimization technique for the weights W is presented, based on proximal descent [START_REF] Parikh | Proximal algorithms[END_REF], and its implementation in a subspace clustering algorithm is then detailed.

Proximal Descent

Studied as a function of variable W only, the cost function is of the general form J(W ) = F (W )+γG(W ), where F corresponds to a differentiable model for subspace clustering and G enforces some constraints on the solutions. This form of functions has recently gained interest in the machine learning community, and proximal descent has thus been studied as an alternative to standard, gradientbased optimization techniques [START_REF] Bach | Optimization with sparsityinducing penalties[END_REF].

As F is differentiable, usual optimization techniques would suggest to use an iterative algorithm based on its gradient, such as gradient descent. This technique considers an update equation of the general form

W t+1 = W t -η • ∇F (W t )
where t is the iteration index and η a descent step size. This simple optimization scheme provides an iterative algorithm in order to minimize any convex function F , starting from any W 0 and iterating until convergence.

Function G is not differentiable, therefore its gradient ∇G is not defined for each W t . Proximal descent [START_REF] Parikh | Proximal algorithms[END_REF] suggests to split the optimization of the two functions F and G, and to enrich gradient descent in the following way:

W t+1 = prox γ L G W t - 1 L ∇F (W t ) (5)
where where L > 0 is a descent step size similar to η. That is, in order to solve a global minimization problem, proximal descent solves a minimization problem as defined by Equation ( 6) at each step of the iteration. The justification for this additional minimization problem is the following: at each step, starting from the current step of iteration W t , proximal descent moves according to -∇F as would any descent approach, then according to prox G . Figure 2 illustrates this principle, showing the space of admissible solutions in blue. This justifies the definition of the proximal operator as well: given an intermediate descent step W t , the solution W to Equation ( 6) is sought for in the neighborhood of W t , hence the first term W -W 2 .

prox γ L G (W ) = argmin W 1 2 W -W 2 + γ L G(W ) (6)
The definition of Equation ( 6) can seem counter-intuitive, as it proposes to solve an additional minimization problem at each step of the algorithm, in order to minimize the global function J. The key ingredient in order to efficiently solve this problem is the notion of proximal operator: a closed-form expression or algorithm to efficiently solve the optimization problem defined by Equation (6).

A Proximal Operator for G

We establish in the following theorem a proximal operator for the penalty term G(W ). Let K denote the vector (1, 1, . . .

1) ∈ R 1×d , such that K • K = d. Theorem 1. Let G r (W r ) = | d p=1 (w rp ) -1| and L ∈ R. prox γ L Gr (W r ) = W r + 1 d K • 1 + prox γd L |•| (K • W r -1) -K • W r (7)
where prox λ|•| (x) = sign(x) max(|x| -λ, 0).

Moreover, prox γ L G (W ) = prox γ L Gr (W r ) r=1...c ∈ R d×c .
Proof. The proof uses results from [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] and [START_REF] Parikh | Proximal algorithms[END_REF]. First, G r (W r ) = φ(K • W r ) where φ(x) = |x-1|. Using the translation and semi-orthogonal linear transform properties [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]:

prox Gr (W r ) = W r + 1 d K • prox φ (K • W r ) -K • W r ) = W r + 1 d K • 1 + prox d|•| (K • W r -1) -K • W r )
Hence the expression of prox γ L Gr by the postcomposition property [START_REF] Parikh | Proximal algorithms[END_REF]. Finally, prox γ L G is computed using the separable sum property of proximal operators [START_REF] Parikh | Proximal algorithms[END_REF].

Equation [START_REF] Wang | K-subspace clustering[END_REF] gives the expression of the proximal operator of the G function. This operator can be used to efficiently implement the scheme defined in Equation ( 5) to find the solution W minimizing function J.

As η for gradient descent, the choice of constant L determines the speed of convergence of the descent. We experimentally observe that setting L to trace(H -1 ) yields good results, where H is the Hessian matrix of F (as a function of W ). As we set v = 2 in Equation ( 2), F is a sum of quadratic functions and its Hessian does not depend on W . It can thus be computed once and for all before running proximal descent.

A Fuzzy Subspace Algorithm: PFSCM

The previous mathematical results allow us to efficiently implement proximal descent for the identification of the subspaces. Combining the expression of prox γ L G with the results of Section 3.2, we introduce the PFSCM algorithm for fuzzy subspace clustering, detailed in Algorithm 1 and commented on below. PFSCM is a FCM-style alternate optimization algorithm which substitutes the exact update of the weights W with proximal descent.

Initialization is a typical issue of k-means-like algorithms. In this paper, initial centers are randomly chosen and each cluster receives uniform weights for all dimensions. As C and U are optimized before W , this is equivalent to running the FCM algorithm and initializing the subspace algorithm with the result, as is often done for subspace clustering algorithms such as Keller and Klawonn [START_REF] Keller | Fuzzy clustering with weighting of data variables[END_REF] or Gustafson-Kessel [START_REF] Gustafson | Fuzzy clustering with a fuzzy covariance matrix[END_REF] algorithms.

As most partitioning algorithms, the number c of clusters to identify must be set by the user, as well as the constant γ > 0. PFSCM is not sensitive to precise values of γ and just needs a value "large enough" (e.g. 10 4 ) in order to impose ∀r, The convergence criteria are defined as the distance between the current and the previous values of the parameters being optimized. In particular, convergence for (C, U ) is defined as C t -C t+1 2 < ∧ U t -U t+1 2 < . PFSCM outputs U , C and . In order to exploit the result of the algorithm, it may be of interest to extract the dimension associated to each cluster. To that aim we propose to post-process the matrix W using a threshold parameter cut to cut out the irrelevant dimensions in a simple fashion: a dimension p for a cluster C r is considered relevant if w rp > cut.

However, instead of using a threshold parameter, an improvement of the algorithm would be to identify the irrelevant dimensions during the descent in order to get rid of them and not have them influence the identification of the centers. This motivates the research of sparsity-inducing models, which are discussed in the next section.

Towards a Sparsity-Inducing Penalty Function

This section discusses extensions of PFSCM aiming at producing sparse solutions: as can be experimentally observed (see Section 6), PFSCM leads to solutions where the weights can be small but different from 0. In order to produce sparse estimations of the subspaces, the previous function J must be changed. In this section, various solutions are discussed to define variants for G, listing some of the difficulties or the similarities with other approaches of the literature.

Indeed, for function G defined in Equation ( 2) and for large values of γ, prox γG acts as a projection towards the simplex ∆ = W r ∈ R d d p=1 w rp = 1 for each W r . The framework presented in this paper would encourage looking for a function which would favor sparse projections on ∆.

1 -norm regularization. The most common sparsity-inducing norm, the 1 norm, is useless in the current setting: indeed, after projection using prox G , the solutions of interest all verify ∀r, W r 1 = 1, and thus this norm does not discriminate any solution under this constraint. A better penalty-inducing function H must be found, which still needs to be "proximable". Moreover, even if G and H are proximable, this is not necessarily the case of their sum G + H [START_REF] Parikh | Proximal algorithms[END_REF].

Negentropy based regularization. An alternative to 1 and other non-differentiable norms may come from negentropy-based regularization, already used in conjunction of standard FCM clustering to regularize the number of clusters (see e.g. [START_REF] Sahbi | Validity of fuzzy clustering using entropy regularization[END_REF]). The study conducted in Section 6 has not yet be extended to negentropy-regularized subspace algorithms. However, Jing et al. [START_REF] Jing | An entropy weighting k-means algorithm for subspace clustering of high-dimensional sparse data[END_REF] use a very similar technique for a slightly different subspace clustering model, in a different context: as they consider sparse data, they use entropy (not negentropy) in order not to promote sparsity, but to promote subspaces of high dimensionality.

Projection-based regularization. Closer to the spirit of the present paper, another direction of work is to reformulate G as an actual projection on ∆ and then to introduce a sparsity-promoting penalty. First define G 1 to be the following penalty function:

G 1 (W r ) = 0 if d p=1 w rp = 1 ∞ otherwise
Function G 1 is convex and thus has a proximal operator:

prox G1 (W r ) = argmin W r ∈R d G 1 (W r ) + 1 2 W r -W r 2 (8)
= argmin

W r ∈∆ 1 2 W r -W r 2 (9) = proj ∆ (W r ) (10)
that is, the proximal operator of G 1 is the projection on ∆. Unlike the function G in Equation ( 2), this projection does not depend on the value of γ: G 1 strictly enforces the constraint, but still does not discriminate between the points of ∆. The goal now becomes to assign a price to the solutions on ∆, in order to favor sparse ones. However, this direction is left for future works.

Gray cluster

Dark cluster Red cluster Blue cluster w 1 w 2 w 1 w 2 w 1 w 2 w 1 w 2 Weights 0.528 0.472 0.063 0.937 0.027 0.973 0.964 0.036 Table 1: Computed weights for the example given in Figure 3. Column w 1 (resp. w 2 ) denotes the weight associated to the x-axis (resp. y-axis).

Experimental Study

The proposed PFSCM algorithm has been tested on artificial data in order to study its ability to correctly identify centers of non-circular clusters, as well as the dimensions that are relevant to describe the clusters. The results show the effectiveness of PFSCM in detecting the clusters and their subspaces. Moreover, PFSCM is compared to Keller and Klawonn's algorithm [START_REF] Keller | Fuzzy clustering with weighting of data variables[END_REF] and shows to provide a better estimation of the dimensionality of the subspaces.

Illustrative Example

This subsection presents illustrative experiments of subspace clustering using the PFSCM algorithm, similar to the visual examples given in [START_REF] Keller | Fuzzy clustering with weighting of data variables[END_REF]. The first example in d = 2 dimensions is represented in Figure 3: four clusters are generated, one of them (the top gray one in Figure 3) being circular while the others have a very low variance in one dimension. PFSCM is run with c = 4, m = 2, and γ = 1000.

In Figure 3 the points are colored according to the cluster C r for which u ri is maximum and Table 1 presents the weights computed for each dimension and cluster. It can be observed that PFSCM correctly identifies the desired clusters and their dimensions: the two weights (w 1 , w 2 ) found for the circular cluster are w 2 ). It is worth noting that, for this specific instance, some points close to the blue cluster are assigned to one of the horizontal clusters, as it minimizes the cost function. This kind of inliers is frequent in subspace clustering problems, and naturally leads to the use of fuzzy membership values (u ri ). It was also studied for its own sake by Guillon et al. [START_REF] Guillon | Laplacian regularization for fuzzy subspace clustering[END_REF], which suggest the use of a Laplacian-based regularization term to prevent it.

Another example is given in Figure 4. Again, 4 clusters are depicted and colored according to the cluster maximizing u ri . These clusters live in various subspaces of various dimensionalities, but are still identified by PFSCM.

Artificial Data Set with Ellipsoidal Clusters

The previous example is generalized to higher dimensions and several algorithms of fuzzy subspace clustering are compared: PFSCM, Keller and Klawonn (KK), Gustafson-Kessel (GK) and WLFC algorithm from [START_REF] Guillon | Laplacian regularization for fuzzy subspace clustering[END_REF].

Experimental Protocol

Considered Data. In order to validate PFSCM, the previous experiment is generalized to higher dimensions, more precisely to artificial data of dimension d ∈ {5, 7, 9, 11, 13, 15}. For each experiment, k = 4 centers are randomly generated in the hypercube [-3, 3] d . Then, d r dimensions j 1 , . . . j dr are randomly picked, with d r randomly chosen between 1 and d -3. Dimensions j 1 , • • • , j dr are thereafter called the "relevant dimensions" for cluster C r . For each cluster, 100 points are generated according to a Gaussian distribution, with variance v j < 0.1 for dimensions j 1 , . . . , j dr and v j ∈ [0.5, 0.9] for other dimensions. The generated points in cluster r in dimension j thus follow X r ∼ N (c r , v j ).

Algorithm Parameters. All algorithms are initialized with FCM centers and use m = v = 2 and c = 4. PFSCM is ran with γ = 1000 and c = 4. All algorithms were implementend using the same convergence criterion, with = 10 -4 . The parameter cut is set to 1 2d , which is a simple rule of thumb to identify the dimensions selected as relevant by the algorithms in each considered dimension d.

Quality Criteria. All algorithms are evaluated on three metrics in order to qualify their results and their ability to discover the desired clusters and subspaces, and their dimensions.

First, let δ = 4 r=1 ĉr -c r 2 be the sum of the Euclidean distances between the generated centers and the computed ones (c r ): this metric is a standard quality criterion for evaluating the produced clusters. A low value means that the computed centers are close to the original ones.

We also consider θ defined as the percentage of clusters for which all relevant dimensions are correctly identified by the algorithm: the relevant dimensions are correctly identified if

w rj > cut ⇔ j ∈ {j 1 , • • • j dr }.
Finally, for the clusters for which the relevant dimensions have been correctly identified, let the weight ratio be φ = ω1 ωj dr where ω 1 is the largest computed weight and ω j dr the smallest computed weight for the relevant dimensions. This metric measures the distortion of the cluster between the relevant dimensions, as estimated by the algorithms. Cluster identification and stability. The results of the experiment are gathered from [START_REF] Guillon | Laplacian regularization for fuzzy subspace clustering[END_REF] and reported in three tables. Table 2 and Table 3 give the means and standard deviations for δ and φ. Table 4 gives the criterion θ, computed over 100 runs of each algorithm. The low mean distances δ in Table 2 show that WLFC and PFSCM succeed in identifying the centers of the generated clusters, while K&K is not as efficient. The GK algorithm performs well in low dimensions, but significantly worse when d increases. Moreover, WLFC and PFSCM appear more stable than K&K and GK, judging by the lower standard deviation.

Dimensionality. The values of θ in Table 4 show that WLFC and K&K succeed in estimating the dimensionality of roughly half of the generated clusters whereas PFSCM performs better when d increases. These results may be related to the choice of the parameter cut used to defuzzify W . This table also shows that GK has troubles identifying the relevant dimensions of the generated clusters; as θ gets very low in higher dimensions, the values of φ for these dimensions is not relevant. A possible explanation of these poor results is that GK is not constrained to search for axes-parallel subspaces and converges towards another solution, as there are not enough points to guide the algorithm. This would also explain the results of the algorithm in Table 2. In Table 3, WLFC, K&K and PFSCM show similar results, bearing in mind the very large standard deviations, PFSCM seeming slightly more stable.

Computation time. In terms of computational time, K&K is significantly faster than all other algorithms. WLFC appears to be globally twice as slow as K&K on average regardless of the dimension. GK is two to five times slower than K&K, probably because of the use of expensive matrix operations. Finally, PF-SCM is much slower than the other algorithms: in higher dimensions it is more than ten times slower than K&K. This may be due to a naive implementation of the inner loop and to the choice of the convergence criterion.

Axes-parallel Gustafson-Kessel algorithm. For the aforementioned reasons, the general GK algorithm is not suited to this experiment. As presented in [START_REF] Borgelt | Fuzzy subspace clustering[END_REF], this algorithm can be modified in order to identify only axes-parallel soft subspaces. We thus ran a specific experiment in order to compare PFSCM to axesparallel GK, which shows that both algorithms perform similarily in terms of the δ metric. However, as GK algorithm describes the subspaces in terms of the variance of their dimension, it provides no indication as to if a dimension is actually relevant or not in the description of a cluster. This feature also seems to preclude its use for the main goal of this work, which is the identification of sparse subspaces.

Conclusion and Future Works

In this paper, a new approach to solve the fuzzy subspace clustering problem with a cost function involving a non-differentiable term has been introduced. The cost function is expressed as the sum of two functions, a model and a penalty or regularizer term. Advanced optimization techniques are explored, which replace the standard update equations of fuzzy c-means-like algorithms.

Experiments on synthetic data show the relevance of the proposed approach, that appears to correctly identify all the relevant dimensions and not more, whereas Keller and Klawonn's algorithm tends to underestimate the number of relevant dimensions. This provides more information about the importance of each dimension for the subspaces and clusters.

Future works will aim at generalizing this approach around the same key ideas based on the discussion initiated in the paper: a differentiable function matching the specification of the problem and one or several penalty functions, expressing constraints on the shape of the solution. The introduction of regularization terms for parameters other than W will also be studied. Finally, more efficient descent schemes will be considered, in order to speed up the descent.

Figure 1 :

 1 Figure 1: Two clusters, contained in two different planes: c 1 in (x, z) and c 2 in (x, y).

Figure 2 :

 2 Figure 2: Principle of the two-steps proximal descent.

d p=1 w rp = 1 .Algorithm 1 :

 11 The algorithm then iterates the update of U , C and W , similarly to alternate optimization in k-means-like algorithms. It consists of two alternate inner loops: the regular parameters C and U are optimized separately from W , which requires the special optimization procedure described in the previous subsection. Variables C and U are optimized one last time at the end of the algorithm, in order to guarantee that the result takes the final computed weights into account.Data: X: data matrix Parameters: c,γ: numbers; Variables: C, U, W: arrays; W last : array Initialization: W r ← (1, 1, . . . 1) for each C r ; C ← random centers Output: C, U, W repeat repeat Update U according to Equation (3); Update C according to Equation (4) until convergence(C, U); repeat Update W according to Equation (6) until convergence(W); W last ←-W until convergence(W last ); Update U and C one last time. Proximal fuzzy subspace clustering PFSCM algorithm

Figure 3 :

 3 Figure 3: Clustering example in two dimensions

Figure 4 :

 4 Figure 4: Clustering example in three dimensions

Table 2 :

 2 Average and standard deviation of δ

	d	WLFC	K&K	GK	PFSCM
	5 0.58 ± 0.09 1.18 ± 0.77 0.59 ± 0.10 0.60 ± 0.11
	7 0.72 ± 0.12 1.49 ± 1.03 0.75 ± 1.67 0.73 ± 0.13
	9 0.86 ± 0.14 1.98 ± 1.58 3.36 ± 5.09 0.84 ± 0.13
	11 0.99 ± 0.13 2.02 ± 1.36 9.78 ± 10.4 0.98 ± 0.12
	13 1.12 ± 0.11 1.77 ± 0.91 19.4 ± 13.0 1.11 ± 0.10

Table 3 :

 3 Average and standard deviation of φ

	d	WLFC	K&K	GK	PFSCM
	5 1.83 ± 1.63 1.59 ± 1.26 11.7 ± 16.7 1.32 ± 1.12
	7 1.93 ± 2.09 1.74 ± 1.83 13.8 ± 18.7 1.73 ± 1.43
	9 2.20 ± 2.64 2.14 ± 2.54 15.2 ± 16.9 2.37 ± 2.10
	11 2.06 ± 2.42 1.97 ± 2.01	NA	1.93 ± 2.04
	13 2.99 ± 3.82 2.99 ± 4.00	NA	3.12 ± 2.19

Table 4 :

 4 Average and standard deviation of θ

	d WLFC K&K GK PFSCM
	5	52%	51% 29%	63%
	7	53%	51% 17%	68%
	9	48%	46% 12%	74%
	11	48%	50%	3%	78%
	13	58%	56%	1%	85%
	6.2.2. Obtained Results

Acknowledgments

The authors would like to thank Nikhil R. Pal for his valuable help and comments regarding this work.

This work was supported by the French "Investissement d'avenir" program (Big-data -Cloud Computing) and was carried out as part of the REQUEST project.


domain: sdv sdv

text:
Isabelle Gardin
email: isabelle.gardin@chb.unicancer.fr

V Grégoire

D Gibon

H Kirisli

D Pasquier

J Thariat

P Vera

Radiomics: Principles and Radiotherapy Applications

Keywords: Radiomics, Radiotherapy, Chemoradiotherapy, Medical Imaging

Radiomics is defined as the extraction of a large quantity of quantitative image features. The different radiomic indexes that have been proposed in the literature are described as well as the various factors that have an impact on the robustness of these indexes. We will see that several hundred quantitative features can be extracted per lesion and imaging modality. The ever-growing number of features studied raises the question of the statistical method of analysis used.

This review addresses the research supporting the clinical use of radiomics in oncology in the staging of disease, discrimination between healthy and pathological tissues, the identification of genetic features, the prediction of patient survival, the response to treatment, the recurrence after radiotherapy and chemoradiotherapy and the side effects.

Based on the existing literature, it remains difficult to identify features that should be used for current clinical practice.

Concept of radiomics

Medical imaging, such as computed tomography (CT), positron emission tomography (PET) using FDG radioactive glucose analogue, and magnetic resonance imaging (MRI) are used routinely in the clinical management of cancer patients. They tend to play a leading role in personalised treatment based on tumour imaging. Personalised medicine has been largely developed using invasive techniques based on genomics and proteomics. However, there is a spatial and temporal heterogeneity of tumour features [START_REF] Marusyk | Intra-tumour heterogeneity: a looking glass for cancer?[END_REF] that makes it difficult to repeat biopsies, both in space and time, where imaging can readily give a spatial representation of the lesion and has the advantage of being easy to repeat.

Radiomics is defined as the extraction of a large quantity of quantitative features from the image [START_REF] Kumar | Radiomics: The process and the challenges[END_REF]. It allows highlighting the heterogeneity, texture, and shape of the lesion corresponding to tissue characteristics such as necrosis, metabolic activity, tissue density, and vascularity [START_REF] Alic | Quantification of heterogeneity as a biomarker in tumor imaging: A systematic review[END_REF]. Medical interpretation of images is generally based on a simple visual interpretation of the contrast. Although this type of interpretation has proven to be very effective in the management of cancer patients, it remains qualitative and subjective, while quantitative features are needed for reasons of objectivity and reproducibility.

There are two ideas underlying the concept of radiomics. The first is that tumour features at the tissue, cellular, and/or genomic level would also be highlighted in medical imaging [START_REF] Yip | Applications and limitations of radiomics[END_REF]. It is considered that image features are strongly correlated with clinical and/or biological features. The second rationale is that the information obtained by the image would be complementary to the other data sources, thus increasing the number of features [START_REF] Lambin | Radiomics: Extracting more information from medical images using advanced feature analysis[END_REF] and the knowledge of the tumour.

The aim is to highlight recent results of the contribution of radiomics to radiotherapy in recent years and to assess the current state of the literature in this field. In the first part of this paper, we will describe the different radiomic features that have been proposed in the literature, the various factors that have an impact on the robustness of these indexes, and the most salient clinical results. We will see that several hundred quantitative features can be extracted by lesion and imaging modality, requiring the use of statistical tools to determine which of these features are relevant. We will conclude with the current limitations of radiomics and the perspectives that this new concept reveals.

Image features

Numerous image features have been proposed in the literature [START_REF] Sollini | PET Radiomics in NSCLC: State of the art and a proposal for harmonization of methodology[END_REF]. Features based on the shape and size of the lesion, histogram features based on first-order statistics, texture features, and filter-and model-based features are generally distinguished (see Figure 1).

Shape-and size-based features

In oncology, it is generally considered that a small tumour volume (TV) leads to a better prognosis than a large lesion. The reduction in lesion volume in CT imaging, assessed by measuring the longest length of the lesion, is a predictive factor for treatment response (RECIST 1.1 criteria [START_REF] Eisenhauer | New response evaluation criteria in solid tumours: revised RECIST guideline (version 1.1)[END_REF]). A review of the literature from Van de Wiele et al. [START_REF] Van De Wiele | Predictive and prognostic value of metabolic tumour volume and total lesion glycolysis in solid tumours[END_REF] lists numerous publications showing the predictive and prognostic value of the TV determined using FDG PET images for solid tumours. Tumour volume is thus an imaging feature of the first choice.

Similarly, in oncology, spherical lesions with sharp edges are frequently considered less aggressive than diffuse and infiltrating tumours. Thus, imaging features corresponding to a metric of this phenomenon have been investigated, such as the sphericity and the extent.

Hofheinz et al. showed the prognostic value of the sphericity on pre-therapeutic FDG-positive lesions in PET imaging in a cohort of 37 patients with head and neck tumours (H&N) [START_REF] Hofheinz | Increased evidence for the prognostic value of primary tumor asphericity in pretherapeutic FDG PET for risk stratification in patients with head and neck cancer[END_REF].

El Naqa et al. showed that the extent, in combination with other FDG PET imaging features, was involved in predicting the response to chemoradiotherapy (CRT) of H&N tumours [START_REF] El Naqa | Exploring feature-based approaches in PET images for predicting cancer treatment outcomes[END_REF].

Histogram-based features

In imaging, the basic feature is the intensity (I) of the voxel. In CT, this corresponds to the Hounsfield units (HU) linked to the tissue density travelled by the X-rays, which is widely used for the delineation of the gross tumour volume (GTV) and the organs at risk (OAR). In PET imaging, I corresponds to a normalised FDG uptake, called SUV for standard uptake value, representing the consumption of glucose by cells based on the fact that tumour cells need more energy than healthy cells. The maximum uptake value, SUVmax, is the most common feature used in PET imaging. FDG-positive tissues are used to delineate the biological target volume proposed for dose painting in RT in the setting of clinical trial [START_REF] Differding | PET imaging biomarkers in head and neck cancer[END_REF].

Furthermore, high FDG uptake areas seen on pre-treatment PET/CT scans can identify tumour subvolumes at greater risk of relapse in patients with non-small cell lung cancer (NSCLC) treated using concomitant CRT [START_REF] Calais | Areas of high 18F-FDG uptake on preradiotherapy PET/CT identify preferential sites of local relapse after chemoradiotherapy for non-small cell lung cancer[END_REF]. In MRI, the interpretation of I is more complex due to the diversity of possible sequences and the multiple signal origins (proton density, relaxation times, diffusion coefficient, vascularity, etc.). In particular, in diffusion-weighted imaging (DWI), the intensity corresponds to the apparent diffusion coefficient (ADC) based on the fact that tumour tissues are composed of anarchic cell proliferation, where water diffusion is usually reduced. On the other hand, dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) allows the visualisation of tumour vascularisation [START_REF] Kumar | Radiomics: The process and the challenges[END_REF][START_REF] Asselin | Quantifying heterogeneity in human tumours using MRI and PET[END_REF].

Descriptive statistics, also called first-order statistics, describe voxel intensities in the volume of interest, here the TV, corresponding to the measure of central tendency and variability. The main imaging feature used to measure the central tendency is the mean (Imean), while the measure of variability includes the standard deviation (ST) or the variance, the coefficient of variation (COV) defined as ST divided by Imean. Bundschuh et al. showed that the tumour heterogeneity of FDG uptake measured by the coefficient of variation was predictive of the response to CRT in rectal cancer [START_REF] Bundschuh | Textural Parameters of Tumor Heterogeneity in 18F-FDG PET/CT for Therapy Response Assessment and Prognosis in Patients with Locally Advanced Rectal Cancer[END_REF].

The frequency histogram is a simple and concise summary of first-order statistical information. Its shape reveals the image homogeneity in the TV. The histogram of a homogeneous lesion is distributed on a narrow peak, whereas a heterogeneous image corresponds to a broad peak of intensities. The general features of intensity distribution are asymmetry (skewness), flattening (kurtosis), energy, and entropy. By analogy with dosevolume histograms used in RT, El Naqa et al. have proposed an intensity-volume histogram to characterise tumour heterogeneity, enabling several indexes to be extracted, such as Ix (the minimum intensity to the x% highest intensity volume) and Vx (the percentage volume with at least x% intensity value) [START_REF] El Naqa | Exploring feature-based approaches in PET images for predicting cancer treatment outcomes[END_REF]. In addition to the feature of the form corresponding to extent, El Naqa et al. showed in their small cohort of patients [START_REF] Hofheinz | Increased evidence for the prognostic value of primary tumor asphericity in pretherapeutic FDG PET for risk stratification in patients with head and neck cancer[END_REF] with H&N cancer treated by CRT that V90 measured on the pre-treatment FDG PET image had the highest univariate predictive power of survival. King et al. [START_REF] King | Head and neck squamous cell carcinoma: Diagnotic performance of diffusion-weighted MR imaging for the prediction of treatment response[END_REF] found in H&N cancer that tumours that responded poorly to CRT showed a significant increase in the mean ADC values on DWI during treatment.

Texture features

The benefit of first-order indexes is that they are simple to calculate when the TV is delineated. However, they do not provide any information about the voxel intensity relative to each other. For this reason, it was proposed to use higher-order statistical methods to carry out a more in-depth analysis of the relationships between voxels, such as texture features. Four main texture matrices have been proposed in the literature. The grey-level co-occurrence matrices (GLCM) characterise the intensity relationships between pairs of neighbouring voxels in all spatial directions [START_REF] Haralick | Textural features for image classification[END_REF]. The grey-level difference matrix (GLDM) [START_REF] Amadasun | Textural features corresponding to textural properties[END_REF] characterises the intensity differences between neighbours. The grey-level run length (GLRLM) [START_REF] Galloway | Texture analysis using gray level run lengths[END_REF] and the grey-level size zone (GLSZM) matrices [START_REF] Thibault | Texture Indexes and Gray Level Size Zone Matrix Application to Cell Nuclei Classification[END_REF] characterise the range of intensities in one direction or in all directions, respectively. The GLRLM features are highly correlated with the GLSZM features and therefore do not provide complementary information [START_REF] Tixier | Intratumor Heterogeneity Characterized by Textural Features on Baseline 18F-FDG PET Images Predicts Response to Concomitant Radiochemotherapy in Esophageal[END_REF]. Thus, in 3D imaging, the GLRL matrix is not of particular interest and is generally not used.

Tixier et al. [START_REF] Tixier | Intratumor Heterogeneity Characterized by Textural Features on Baseline 18F-FDG PET Images Predicts Response to Concomitant Radiochemotherapy in Esophageal[END_REF] studied the predictive value of the response to CRT of 38 texture features in addition to first-order features in a cohort of 41 patients with oesophageal cancer. Their results showed that features of textures enable to differentiate patients better than those of the first order.

Filter-or model-based features

Other approaches based on prior image transformation have been proposed in the literature, mainly on CT and MRI data, defined in spatial domains such as filters (Gabor filters, Law's filters, and Laplacian transform of Gaussian filters (LoG filters) or in the frequency domains (Fourier and wavelet transform) [START_REF] Kumar | Radiomics: The process and the challenges[END_REF]. Other features derived from complex mathematical models have also been proposed in the literature, such as fractal dimension measuring image complexity [START_REF] Lopes | Fractal and multifractal analysis: A review[END_REF].

In order to study the predictive value of distant metastasis (DM) and the prognostic values of radiomic features, Coraller et al. [START_REF] Coroller | CT-based radiomic signature predicts distant metastasis in lung adenocarcinoma[END_REF] extracted 635 indexes from planning CT in 98 patients treated using CRT for lung cancer for analysis and 84 patients for validation. Of the 635 features, 180 were extracted from LoG filters and 384 from wavelet transform maps. The authors found 12 prognostic features and 35 predictive features of DM occurrence, including several LoG filters and wavelet transform features.

Medical applications

Various reviews of the literature [START_REF] Marusyk | Intra-tumour heterogeneity: a looking glass for cancer?[END_REF][START_REF] Kumar | Radiomics: The process and the challenges[END_REF][START_REF] Yip | Applications and limitations of radiomics[END_REF][START_REF] Sollini | PET Radiomics in NSCLC: State of the art and a proposal for harmonization of methodology[END_REF][START_REF] Scalco | Texture analysis of medical images for radiotherapy applications[END_REF] show that many studies support the clinical interest of radiomics in oncology, both in MR, CT, and FDG PET imaging, particularly for patients treated using radiotherapy for many localisations (NSCLC, H&N, oesophagus, breast, cervix, etc.). Before treatment, this concerns the staging of the disease and the identification of genetic features. For treatment planning, this concerns the tissue discrimination between healthy and pathological tissues, the delineation of the GTV, or a boost. For therapeutic monitoring, this concerns the prediction of the patient's survival, the response to treatment, and the occurrence of distant metastases and side effects.

Image features can be extracted from diagnostic or patient monitoring images, as well as from CT [START_REF] Coroller | CT-based radiomic signature predicts distant metastasis in lung adenocarcinoma[END_REF] or MR treatment planning images [START_REF] Shiradkar | Radiomics based targeted radiotherapy planning (Rad-TRaP): A computational framework for prostate cancer treatment planning with MRI[END_REF]. Fave et al. [START_REF] Fave | Can radiomics features be reproducibly measured from CBCT images for patients with non-small cell lung cancer?[END_REF] also showed in 10 NSCLC patients treated using RT that a radiomic analysis could be performed from cone-beam CT (CBCT) images using medical Linac on-board imaging systems. Of 68 features studied, 37 were robust to noise and poor image quality.

Genetic features

Radiomic analysis was used to describe different tumour phenotypes corresponding to distinct genotypes. For instance, Buch et al. [START_REF] Buch | Using Texture Analysis to Determine Human Papillomavirus Status of Oropharyngeal Squamous Cell Carcinomas on CT[END_REF] showed, on contrast-enhanced CT, statistically significant differences in some texture features between human papillomavirus-positive and human papillomavirus-negative oropharyngeal tumours.

Segmentation

Nailon et al. [START_REF] Nailon | Characterisation of radiotherapy planning volumes using textural analysis[END_REF] have proposed an automatic segmentation of the GTV and OAR in bladder cancer based on 27 features of CT imaging. The authors showed that this approach offers significant accuracy using only 3 texture features. Yu et al. combined a radiomic analysis of PET and CT images in H&N tumours, highlighting that some image features had a better discriminating power between healthy and pathological tissues than experts, the UH values (CT), or the SUV (PET imaging) [START_REF] Yu | Coregistered FDG PET/CT-based textural characterization of head and neck cancer for radiation treatment planning[END_REF][START_REF] Yu | Automated Radiation Targeting in Head-and-Neck Cancer Using Region-Based Texture Analysis of PET and CT Images[END_REF].

Radiomic features can also be used for boost planning. Shiradkar et al. [START_REF] Shiradkar | Radiomics based targeted radiotherapy planning (Rad-TRaP): A computational framework for prostate cancer treatment planning with MRI[END_REF] developed a radiomic-based software platform for the treatment of prostate cancer using brachytherapy and external beam radiation therapy (Rad-TRaP, radiation treatment planning for prostate). This platform allows the extraction of 154 radiomic features from T2-weighted and ADC MR images. From a database of 11 patients, the authors identified 11 features (5 from the T2weighted image and 6 from the ADC image), making it possible to discriminate, within the gland, the voxels belonging to the lesion and to calculate the probability of the presence of the tumour for each voxel to define a boost from this probability map.

Prediction of treatment response and recurrence of disease

Aerts et al. [START_REF] Aerts | Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach[END_REF] studied the prognostic value of 440 radiomic features (first-order, form, and texture features (GLCM, GLRLM, and wavelets)) extracted from CT images on 3 cohorts of patients corresponding to a total of 1019 patients with lung cancer or H&N tumours. The authors showed that the results obtained using a radiomic signature generated from 4 features of the image were better than those obtained using TNM staging or the TV alone. The 4 features were energy (first-order statistics), compactness (shape feature), and 2 features of grey-level non-uniformity (texture and wavelets, respectively). In addition, this signature appeared to be associated with gene expression profiles.

In a database of 92 patients, Li et al. [START_REF] Li | Imaging features from pretreatment CT scans are associated with clinical outcomes in nonsmall-cell lung cancer patients treated with stereotactic body radiotherapy[END_REF] showed that features extracted from pre-treatment planning CT have a prognostic value of recurrence after stereotactic body radiotherapy (SBRT) for NSCLC. Because of the difficulty of differentiating between recurrence of lung cancer after stereotactic radiotherapy treatment and benign radiation-induced lung injury (RILI), Mattonen et al. [START_REF] Mattonen | Detection of Local Cancer Recurrence after Stereotactic Ablative Radiation Therapy for Lung Cancer: Physician Performance Versus Radiomic Assessment[END_REF] showed the contribution of radiomics in the rapid prediction of the recurrence of cancerous disease.

Gnep et al. [START_REF] Gnep | Haralick textural features on T2-weighted MRI are associated with biochemical recurrence following radiotherapy for peripheral zone prostate cancer[END_REF] showed that pre-therapeutic features extracted from the co-occurrence matrix of T2-weighted images in patients treated with RT for prostate cancer were strongly correlated with a biochemical recurrence of the disease.

Prediction of side effects

Radiomic analysis has also been proposed to predict the occurrence of side effects due to ionising radiation on OAR in NSCLC and oesophageal cancers, as well as salivary glands and parotids in the treatment of H&N tumours using RT [START_REF] Scalco | Texture analysis of medical images for radiotherapy applications[END_REF].

Moran et al. [START_REF] Moran | Radiomics-based Assessment of Radiation-induced Lung Injury After Stereotactic Body Radiotherapy[END_REF] showed that CT-based radiomic features were significantly correlated with radiation oncologist-scored post-SBRT lung injury, suggesting the potential of radiomics to measure post-SBRT lung injury. Cunliffe et al. [START_REF] Cunliffe | Lung texture in serial thoracic computed tomography scans: correlation of radiomics-based features with radiation therapy dose and radiation pneumonitis development[END_REF] showed, in a database of 106 patients treated with RT for their oesophageal cancer, a significant change in 20 radiomic features between the pre-RT CT scan and the post-RT CT scan as a function of the radiation dose.

Twelve features changed significantly in patients who developed radiation pneumonitis (RP), allowing the authors to develop a classifier based on these CT features. The same team [START_REF] Anthony | Incorporation of pre-therapy 18F-FDG uptake data with CT texture features into a radiomics model for radiation pneumonitis diagnosis[END_REF] showed the additive value of the SUV standard deviation on pre-RT PET images indicating the potential of a multimodality radiomic analysis. on a prospective study, Jahangiri et al. [START_REF] Jahangiri | A prospective study of the feasibility of FDG-PET/CT imaging to quantify radiation-induced lung inflammation in locally advanced nonsmall cell lung cancer patients receiving proton or photon radiotherapy[END_REF] showed the feasibility of FDG PET/CT imaging to quantify radiation-induced lung inflammation in locally advanced NSCLC receiving proton or photon radiotherapy. This study suggests less induction of inflammatory response in both ipsilateral and contralateral lungs of patients treated with proton compared to photon or combined proton-photon RT.

From a radiomic study conducted on a retrospective treatment planning CT database of 249 patients treated using RT for H&N tumours, van Dijk et al. [START_REF] Van Dijk | CT image biomarkers to improve patientspecific prediction of radiation-induced xerostomia and sticky saliva[END_REF] showed the benefits of radiomic analysis in predicting xerostomia and sticky saliva 12 months after RT. With the same goal, Pota et al. [START_REF] Pota | Early prediction of radiotherapy-induced parotid shrinkage and toxicity based on CT radiomics and fuzzy classification[END_REF] developed a model for classifying patients at risk of developing parotid gland shrinkage and 12-month xerostomia. Features came from several information sources (clinical data, dosimetry, and CT).

Factors influencing radiomic features

The very encouraging results of radiomics, however, raise the issue of the robustness of features with respect to the experimental conditions and image processing performed before their extraction as well as the manner of extraction.

Feature definition

In the literature, there are differences in the definition of radiomic features due to a lack of consensus [START_REF] Orlhac | Tumor Texture Analysis in 18F-FDG PET: Relationships Between Texture Parameters, Histogram Indices, Standardized Uptake Values, Metabolic Volumes, and Total Lesion Glycolysis[END_REF]. Sometimes, 2 texture features with different definitions are referred to using the same name. Features are also found under different names, but with the same mathematical expression. Calculation methods can also be different. Some texture features can be extracted from a single slice of the tumour (2D texture) or from the entire tumour volume (3D texture) [START_REF] Ng | Assessment of tumor heterogeneity by CT texture analysis: Can the largest cross-sectional area be used as an alternative to whole tumor analysis?[END_REF][START_REF] Fave | Preliminary investigation into sources of uncertainty in quantitative imaging features[END_REF]. Similarly, some authors considered the different directions of the GLC matrices individually, while others preferred to extract features from the average matrix in 13 directions [START_REF] Hatt | 18F-FDG PET Uptake Characterization Through Texture Analysis: Investigating the Complementary Nature of Heterogeneity and Functional Tumor Volume in a Multi-Cancer Site Patient Cohort[END_REF]. These findings show that an effort to harmonise and define features is still essential [START_REF] Sollini | PET Radiomics in NSCLC: State of the art and a proposal for harmonization of methodology[END_REF][START_REF] Buvat | Tumor Texture Analysis in PET: Where Do We Stand?[END_REF][START_REF] Lambin | Radiomics: The bridge between medical imaging and personalized medicine[END_REF].

Image acquisition and processing

Several authors have studied the robustness of radiomic features according to test-retest studies [START_REF] Balagurunathan | Reproducibility and Prognosis of Quantitative Features Extracted from CT Images[END_REF][START_REF] Balagurunathan | Test-retest reproducibility analysis of lung CT image features[END_REF][START_REF] Tixier | Reproducibility of Tumor Uptake Heterogeneity Characterization Through Textural Feature Analysis in 18F-FDG PET[END_REF][START_REF] Desseroit | Reliability of PET/CT Shape and Heterogeneity Features in Functional and Morphologic Components of Non-Small Cell Lung Cancer Tumors: A Repeatability Analysis in a Prospective Multicenter Cohort[END_REF], modifying acquisition [START_REF] Fave | Preliminary investigation into sources of uncertainty in quantitative imaging features[END_REF], image reconstruction parameters [START_REF] Yan | Impact of Image Reconstruction Settings on Texture Features in 18F-FDG PET[END_REF][START_REF] Galavis | Variability of textural features in FDG PET images due to different acquisition modes and reconstruction parameters[END_REF], and lesion segmentation [START_REF] Hatt | Prognostic value of 18F-FDG PET image-based parameters in oesophageal cancer and impact of tumour delineation methodology[END_REF]. These studies showed that some features are less robust than others, leading some authors to eliminate the least robust features from the radiomic analysis.

These results show that it is necessary to harmonise image acquisition and perform the most automatic image processing and segmentation as possible to increase the robustness and relevance of the radiomic analysis [START_REF] Larue | Quantitative radiomics studies for tissue characterization: A review of technology and methodological procedures[END_REF].

Lesion size

Brooks et al. [START_REF] Brooks | The Effect of Small Tumor Volumes on Studies of Intratumoral Heterogeneity of Tracer Uptake[END_REF] studied the variation in entropy (GLCM) as a function of the number of voxels in TV based on a PET database of 70 patients with cervical cancer. The results showed that entropy is strongly correlated with TV if the latter is less than 45 cm 3 (700 voxels).

Below this threshold, the index was 5 times more sensitive to volume differences. The authors have suggested that it is not appropriate to study this feature for a tumour with a volume less than 45 cm 3 , which in oncology corresponds to a large volume and eliminates a significant number of patients.

Hatt et al. [START_REF] Hatt | 18F-FDG PET Uptake Characterization Through Texture Analysis: Investigating the Complementary Nature of Heterogeneity and Functional Tumor Volume in a Multi-Cancer Site Patient Cohort[END_REF] studied the influence of tumour volume on 4 texture features using PET images based on a database of 555 patients with different cancers (lung, breast, and others).

The authors found that the correlation between the radiomic features and the TV tends to decrease when TV increases.

Resampling image intensity values

An important methodological factor is the image discretisation corresponding to resampling image intensity values. When the number of intensity levels in the TV is large, resampling is necessary. Its aim is to group the very large number of intensity levels into a limited number. This step reduces the computation time of the texture matrices, as well as the noise in the image, but to the detriment of a loss of information [START_REF] Tixier | Reproducibility of Tumor Uptake Heterogeneity Characterization Through Textural Feature Analysis in 18F-FDG PET[END_REF].

Two approaches have been proposed in the PET imaging literature. The first is based on relative resampling using a fixed number (e.g., 16, 32, or 64) of discrete resampled values or bins to divide the SUV range into equally spaced intervals before calculating textural features [START_REF] Tixier | Intratumor Heterogeneity Characterized by Textural Features on Baseline 18F-FDG PET Images Predicts Response to Concomitant Radiochemotherapy in Esophageal[END_REF][START_REF] Orlhac | Tumor Texture Analysis in 18F-FDG PET: Relationships Between Texture Parameters, Histogram Indices, Standardized Uptake Values, Metabolic Volumes, and Total Lesion Glycolysis[END_REF]. Orlhac et al. showed that texture features are not very sensitive to the number of intensities used when the intensity level is greater than or equal to 32 [START_REF] Orlhac | Tumor Texture Analysis in 18F-FDG PET: Relationships Between Texture Parameters, Histogram Indices, Standardized Uptake Values, Metabolic Volumes, and Total Lesion Glycolysis[END_REF]. Hatt et al. concluded that this type of resampling by more than 64 intensity levels did not provide any additional information and advised limiting the number of resampling to this level [START_REF] Hatt | 18F-FDG PET Uptake Characterization Through Texture Analysis: Investigating the Complementary Nature of Heterogeneity and Functional Tumor Volume in a Multi-Cancer Site Patient Cohort[END_REF].

More recently, an absolute resampling process was proposed to resample SUV with a fixed bin size in units of SUV (for example, 0.1 and 0.5), maintaining a constant intensity resolution across all tumour images [START_REF] Leijenaar | The effect of SUV discretization in quantitative FDG-PET Radiomics: The need for standardized methodology in tumor texture analysis[END_REF][START_REF] Orlhac | 18F-FDG PET-derived textural indices reflect tissue-specific uptake pattern in non-small cell lung cancer[END_REF]. It has been shown in PET imaging that absolute resampling has a better predictive power than relative resampling [START_REF] Brooks | The Effect of Small Tumor Volumes on Studies of Intratumoral Heterogeneity of Tracer Uptake[END_REF][START_REF] Leijenaar | The effect of SUV discretization in quantitative FDG-PET Radiomics: The need for standardized methodology in tumor texture analysis[END_REF][START_REF] Orlhac | 18F-FDG PET-derived textural indices reflect tissue-specific uptake pattern in non-small cell lung cancer[END_REF][START_REF] Desbordes | Predictive value of initial FDG-PET features for treatment response and survival in esophageal cancer patients treated with chemo-radiation therapy using a random forest classifier[END_REF].

Methods of radiomic analysis

The ever-growing number of features studied raises the question of the statistical method of analysis used and its quality, especially as patient databases can be limited.

Conventional statistical analysis

The majority of the first articles from the literature studying the predictive value of radiomic features were based on the Mann-Whitney U-test when studying the ability of the feature to differentiate patients into 2 classes (for example, responding and not responding to treatment) or the Kruskal-Wallis test for more than 2 classes. When the feature corresponds to a continuous numerical value, the threshold value of the feature used to optimise the classification of the population into groups is determined using receiver operating characteristic (ROC) curve analysis to calculate the area under the curve, the sensitivity, and the specificity. For prognostic studies of patient survival, Kaplan-Meier's univariate analysis combined with a log-rank test followed by Cox's multivariate analysis is generally used.

The number of features extracted from the image does not necessarily mean that the amount of information increases. Many authors have shown that some features are highly correlated using ranking tests and should be excluded from radiomic analysis based on conventional statistical analysis, leading to a reduction in the number of characteristics studied [START_REF] Orlhac | Tumor Texture Analysis in 18F-FDG PET: Relationships Between Texture Parameters, Histogram Indices, Standardized Uptake Values, Metabolic Volumes, and Total Lesion Glycolysis[END_REF][START_REF] Desbordes | Predictive value of initial FDG-PET features for treatment response and survival in esophageal cancer patients treated with chemo-radiation therapy using a random forest classifier[END_REF].

Orlhac et al. [START_REF] Orlhac | 18F-FDG PET-derived textural indices reflect tissue-specific uptake pattern in non-small cell lung cancer[END_REF] showed on PET images that there was a significant correlation between texture features and TV when a relative resampling is used. This correlation does not exist with absolute resampling, but there is a correlation between texture features and SUVmax. To further reduce the number of features, only those that are the most robust with regard to the experimental conditions and image processing are used [START_REF] Balagurunathan | Reproducibility and Prognosis of Quantitative Features Extracted from CT Images[END_REF]. Despite these strategies for limiting the number of features, radiomic studies based on conventional statistics present many methodological pitfalls and limitations when handling a large number of indexes [START_REF] Chalkidou | False discovery rates in PET and CT studies with texture features: A systematic review[END_REF].

Machine learning methods

Given the number of features, the use of machine learning methods is of great interest. Based on a learning database, the goal of machine learning is to develop algorithms that generate a decision rule, also called a model, to predict new patient outcomes. Numerous machine learning algorithms have been proposed in the literature (see Table 1). Examples include logistic regression, multivariate Cox proportional hazards regression models [START_REF] Aerts | Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach[END_REF], methods based on Bayesian statistics, support vector machine searching the optimal hyperplane separating patient outcomes [START_REF] Cortes | Support-Vector Networks[END_REF], as well as the decision tree method and its extension to random forests [START_REF] Breiman | Random forests[END_REF]. Another group of methods corresponds to Artificial Neural Network and deep learning algorithms simulating neuronal functioning [START_REF] Hosny | Deep learning for lung cancer prognostication: A retrospective multi-cohort radiomics study[END_REF].

The advantages of machine learning methods are that they were created, from the start, to handle the processing of a large number of features. Although they require a learning database, which is quite important because the number of features studied is so large, some are better suited than others to the low ratio between the number of patients and the number of features. Depending on the method, they make few or no assumptions about the relationship between feature and prediction, nor about linear or rank relationships between features.

Many feature selection methods have been proposed in the literature to remove features that are highly correlated [START_REF] Orlhac | Tumor Texture Analysis in 18F-FDG PET: Relationships Between Texture Parameters, Histogram Indices, Standardized Uptake Values, Metabolic Volumes, and Total Lesion Glycolysis[END_REF] or that have limited contribution [START_REF] Desbordes | Predictive value of initial FDG-PET features for treatment response and survival in esophageal cancer patients treated with chemo-radiation therapy using a random forest classifier[END_REF][START_REF] Parmar | Machine Learning methods for Quantitative Radiomic Biomarkers[END_REF] and avoid overfitting.

Here, unlike conventional statistical methods, highly correlated features, in the sense of a rank test, can provide additional information.

These methods are recent and widely used in all areas of big data processing. Therefore, they have a strong potential for methodological development in the future. Their diversity allows clinicians to make a choice depending on the problem, such as the balance between the different classes, the presence of features of different formats (binary, discrete, or continuous), and others.

Ypsilantis et al. [START_REF] Ypsilantis | Predicting Response to Neoadjuvant Chemotherapy with PET Imaging Using Convolutional Neural Networks[END_REF] studied the use of several machine learning methods (ANN, SVM, and RF, among others) to predict the response to neoadjuvant chemotherapy in 107 patients with oesophageal cancer. For this purpose, 103 radiomic features were extracted from FDG PET images. ANN gave the best performances for predicting treatment response with a good classification accuracy of 73.4% ± 5%, followed by RF (65.7% ± 6%) and SVM (60.5% ± 8%). The results were superior to those using SUVmax alone, with an accuracy of 41.0% ± 5%.

From a database of 65 patients with oesophageal cancer treated using CRT and 61 clinical and FDG PET features, Desbordes et al. [START_REF] Desbordes | Predictive value of initial FDG-PET features for treatment response and survival in esophageal cancer patients treated with chemo-radiation therapy using a random forest classifier[END_REF] showed the superiority of RF over conventional statistical analysis. The RF classifier improved the predictive and prognostic values compared to the Mann-Whitney U-test and the univariate Kaplan-Meier survival analysis when applied to several tens of features in a limited patient database.

Conclusions and perspectives

Numerous articles support the clinical use of radiomic analysis in oncology, as well as on MRI, CT than on FDG PET images, in staging of the disease, significantly differentiating the early and advanced stages of the disease, tissue discrimination between healthy and pathological tissues, demonstrating genetic characteristics, predicting the response to treatment, patient survival, and side effects. These results from the literature mainly concern solid tumours for many localisations (pulmonary, H&N, oesophagus, breast, and cervix, among others). They potentially concern any type of treatment (chemotherapy, radiotherapy, and concomitant treatments).

This study of the literature reveals very important perspectives within the context of personalised medicine and highlights many scientific issues not yet resolved, which leaves a large margin of progression for the concept of radiomics.

Numerous factors have been reported to influence the robustness of features, such as acquisition conditions, data reconstruction, tumour volume segmentation, and image resampling. Radiomics requires a harmonisation of practises, which raises the issue of the quality of the database, especially multicentric databases, in a context where the results must be robust to be generalised. Some features are more robust than others, but there is currently no real consensus to establish a list of preferred features or tumour signatures to separate lesions of good and poor prognosis.

Given the number of features studied, the size of databases and their constitution remain issues. In this context, conventional statistical analyses seem too limited and many methodological pitfalls and limitations were described [START_REF] Leijenaar | The effect of SUV discretization in quantitative FDG-PET Radiomics: The need for standardized methodology in tumor texture analysis[END_REF] that may explain the lack of the emergence of radiomic signatures. This is likely due to the multiplication of features. It can be noted that very few articles currently address the contribution of multimodality imaging, genomics, and proteomics, corresponding to very important fields of investigation. These issues can be addressed using machine learning methods. Although many methods of selecting features have been proposed in the literature to reduce the number of features, machine learning methods make it possible to manage a large number of features. They can also incorporate correlated features according to a rank test as being relevant information.

Despite the enthusiasm for radiomics, it requires the definition of robust signatures to be used in clinical routine. The ultimate goal of such an approach would be to use these signatures individually for each patient. These radiomic signatures must be coupled with other relevant information sources (clinical, biological, genomics and biophathological) in order to improve the prognostic and predictive value of the treatment towards personalized medicine. Table Table 1: Main machine learning algorithms used in radiomics with their advantages and limitations.

Algorithm Advantages Limitations

Logistic regression [START_REF] Ypsilantis | Predicting Response to Neoadjuvant Chemotherapy with PET Imaging Using Convolutional Neural Networks[END_REF] -easy to interpret -Does not require input features to be scaled

FigureFigure 1 :

 1 Figure Figure 1: Different categories of radiomic features and non-exhaustive list of associated features. Example with two lung lesions in PET/CT imaging. The filter-based example corresponds to wavelet decomposition.

-

  Does not require any tuning -Easy to regularize -Easy to implement and very efficient to train -A good baseline to compare the performances with those of other more complex algorithms -Cannot solve non-linear problems -Features must be independent -Risk of overfitting Multivariate Cox proportional hazards regression model [30] -Does not require the division of the studied features into discrete classes (keep the information as precise as possible) -Limited by its very restrictive application conditions (hazards should be proportional) proper network structure: fine tuning of the model hyperparameters

Financial disclosure statement

This work was supported by grant from the Cancéropôle Nord-Ouest, Farnce