Data augmentation and feature selection for automatic model recommendation in computational physics

Classification algorithms have recently found applications in computational physics for the selection of numerical methods or models adapted to the environment and the state of the physical system. For such classification tasks, labeled training data come from numerical simulations and generally correspond to physical fields discretized on a mesh. Three challenging difficulties arise: the lack of training data, their high dimensionality, and the non-applicability of common data augmentation techniques to physics data. This article introduces two algorithms to address these issues, one for dimensionality reduction via feature selection, and one for data augmentation. These algorithms are combined with a wide variety of classifiers for their evaluation. When combined with a stacking ensemble made of six multilayer perceptrons and a ridge logistic regression, they enable reaching an accuracy of 90% on our classification problem for nonlinear structural mechanics.


Introduction
Classification problems can be encountered in various disciplines such as handwritten text recognition [1], document classification [2], and computer-aided diagnosis in the medical field [3], among many others. In numerical analysis, classification algorithms are getting more and more attention for the selection of efficient numerical models that can predict the behavior of a physical system with very different states or under various configurations of its environment [4,5,6,7,8,9,10,11]. Classifiers have been used as reduced-order model (ROM) selectors in [4,5,8,9,11] in computational mechanics, enabling the computation of approximate solutions at lower cost by replacing a generic high-fidelity numerical model by a specific (or local) ROM adapted to the simulation's context. Reduced-order modeling [12,13] consists in identifying an appropriate low-dimensional subspace on which the governing equations are projected in order to reduce the number of degrees of freedom of the solution. In [11], the combination of a classifier with a dictionary of local ROMs has been termed dictionary-based ROM-net. Such approaches are promising numerical methods using both physics equations and a collection of latent spaces to compute approximations of solutions lying in nonlinear manifolds.
Dictionary-based ROM-nets use a physics-informed automatic data labeling procedure based on the clustering of numerical simulations. Due to the cost of numerical simulations, training examples for classification are limited in number. Moreover, the dimensionality of input data can be very high, especially when dealing with physical fields discretized on a mesh (finitedifference methods [14], finite-element method [15], finite-volume method [16]) or with bond graphs modeling engineering systems [17].
When classification data are high-dimensional, dimensionality reduction techniques can be applied to reduce the amount of information to be analyzed by the classifier. For classification problems where the dimension of the input data is higher than the number of training examples, dimensionality reduction is crucial to avoid overfitting. In addition, when considering physical fields discretized on a mesh, the dimension of the input space can reach 10 6 to 10 8 for industrial problems. In such cases, the input data are too hard to manipulate, which dramatically slows down the training process for the classifier and thus restrains the exploration of the hyperparameters space, as it requires multiple runs of the training process with different values for the hyperparameters. Applying data augmentation techniques to increase the number of examples in the training set is also impossible, as it would cause memory problems. Therefore, dimensionality reduction is recommended not only for reducing the risk of overfitting, but also for facilitating the training phase and enabling data augmentation.
Feature selection [18] aims at decreasing the number of features by selecting a subset of the original features. It differs from feature extraction, where new features are created from the original ones (e.g. Principal Component Analysis, PCA, and more generally encoders taken from undercomplete autoencoders [19]). Feature selection can be seen as applying a mask to a high-dimensional random vector to get a low-dimensional random vector containing the most relevant information. It is preferred over autoencoders when interpretability is important [20]. Furthermore, contrary to undercomplete autoencoders trained with the mean squared error loss, most feature selection algorithms do not intend to find reduced features enabling the reconstruction of the input: features are selected for the purpose of predicting class labels, which makes these algorithms more goal-oriented for supervised learning tasks.
Among the existing feature selection algorithms, univariate filter methods consist in computing a score for each feature and ranking the features according to their scores. The score measures how relevant a feature is for the prediction of the output variable. If N f is the target number of features, then the N f features with the highest scores are selected, and the others are discarded. The major drawback of univariate filter methods is that they do not account for relations between the selected features. The resulting set of selected features may then contain redundant features. To address this issue, the minimum redundancy maximum relevance (mRMR) algorithm [21,22] tries to find a tradeoff between relevance and redundancy. However, for very large numbers of features like in computational physics, evaluating the redundancy is very computationally demanding. Fortunately, working on physics data provides other possibilities to define a redundancy measure. In this paper, we propose a new feature selection algorithm suitable for features coming from the same physical quantity but corresponding to different points in a space-time discretization. It is assumed that this physical quantity, defined as a function of space and/or time, has some smoothness properties. This is often the case in physics, where the physical quantity satisfies partial differential equations and boundary conditions. In [23], it is shown that the solution of Poisson's equation on a Lipschitz domain in R 3 with a L 2 source term and Dirichlet or Neumann boundary conditions is continuous. Poisson's equation is well-known in physics, and can be found for example in electrostatics, in Gauss's law for gravity, in the stationary heat equation, and in the stationary particle diffusion equation. If the features of a random vector contain the discretized values of a smooth function of space and time, then their correlations are related to their proximities on the space-time grid. The approach presented in this paper is depicted as a geostatistical variant of mRMR algorithm, in the sense that it consists in modeling the redundancy as a function of space and time.
Once the dimension of the input space is reduced, another challenge of the classification problems encountered in computational physics must be addressed: the lack of training data. Data augmentation refers to techniques aiming at enlarging the training set by generating new examples from the original ones. For image classification, many class-preserving operations can be used to create new images, such as translations, rotations, cropping, scaling, and changes in colors, brightness and contrast. Unfortunately, these common techniques cannot be used when considering physics data. For this type of data, new examples can be generated using generative adversarial networks (GAN [24], see [25] for the use of deep convolutional GANs in computational fluid dynamics). However, training GANs is quite complex in practice and may also be made more difficult by the lack of training examples. More simply, new data can be generated by convex combinations of the original examples. SMOTE [26] takes convex combinations of input data with their nearest neighbors in the input space. ADASYN [27] uses the same idea but focuses more on examples that are hard to learn, i.e. those having examples of a foreign class in their neighborhoods. Both data augmentation algorithms use k-nearest neighbors algorithm and thus compute Euclidean distances in the input space. When working on high-dimensional physics data, this approach may suffer from the curse of dimensionality [28]. In addition, defining neighborhoods with the Euclidean distance in the input space is not always appropriate, since dictionary-based ROM-nets use physics-aware dissimilarities to label the data, such as distances on the primal variable or on a quantity of interest. The data augmentation algorithm developed in this article consists in growing sets around original examples by incrementally adding nearest neighbors in terms of the dissimilarity measure used for the automatic data labeling procedure. These sets are used to generate new data by convex combinations. Contrary to SMOTE and ADASYN, the risk of generating new data with wrong labels is controlled by checking that the convex hulls of the growing sets do not contain any example belonging to a foreign class.
In sum, the contributions of this paper are motivated by difficulties encountered in our previous work on ROM-nets [11]. These difficulties are inherent to classification tasks on simulation data and can be summarized in three main issues: • the lack of training data due to the expensive data labeling procedure involving simulations with a high-fidelity model (risk of overfitting); • the high dimensionality of input data (risk of overfitting); • most common data augmentation techniques are not applicable to physics data.
The feature selection and data augmentation strategies introduced in this paper are developed to tackle these difficulties. Classification problems encountered in computational physics are described in Section 2. Section 3 presents the classification problem studied in this paper. The feature selection algorithm is described in Section 4 and is shown to efficiently remove irrelevant and redundant features. Section 5 presents the data augmentation algorithm, which successfully generates a large amount of new data with correct labels. Finally, Section 6 evaluates both algorithms in conjunction with 14 different classifiers. On our classification task, the average accuracy gain due to data augmentation is 4.98%. Using ensemble methods on classifiers combined with our algorithms enables reaching a classification accuracy of 90%.
2 Classification in the context of numerical modeling 2.1 Classification: a brief review Supervised learning is the task of learning the correspondence between input data X and outputs Y from a training set of input-output pairs {(x i , y i )} 1≤i≤N . Supervised machine learning problems fall into two categories: regression problems, for which the outputs take continuous values, and classification problems, consisting in the prediction of categorical labels. This paper focuses on the latter, with the additional assumptions that X is a continuous multivariate random variable having a probability density function p X : X → R + , and that any observation x ∈ X is associated to a single label y. The discrete random variable Y follows a categorical distribution (or multinoulli distribution) whose probability mass function is defined by: where K is the number of categories (or classes), δ is the Dirac delta function, and P Y (k) denotes the probability of the event Y = k for a given label k ∈ [[1; K]]. The labeled training data are drawn from the joint probability distribution p X,Y , called the data-generating distribution. As X is continuous and Y is discrete, p X,Y is a mixed joint density and can be obtained with the formula: with p X Y being the class-conditional probability distribution.
In the present paper, we are interested in single-label multiclass problems. Hence, the classification problem considered here reads: given an integer K ≥ 2 and a training set K]] to assign any observation x ∈ X to the correct class, with θ denoting the parameters of the classifier. However, reaching the highest possible accuracy on the training set is not the objective to be pursued, since it usually leads to overfitting. Indeed, the classifier is supposed to be applied to new unseen data, or test data, after the training phase. Therefore, the generalization ability of the classifier is at least as important as its performance on the training set. A classifier with high capacity 1 perfectly fits training data but is very sensitive to noise, leading to high test error and thus overfitting. On the other hand, a classifier with low capacity can produce smaller error gaps between training and test predictions, but such a classifier may not be able to fit the data, which is called underfitting. This dilemma is known as the bias-variance tradeoff : low model capacity leads to high bias, while high model capacity leads to high variance.
For a given observation x ∈ X , probabilistic classification algorithms estimate the membership probabilities P model (y x; θ) for each class y ∈ [[1; K]]. The classifier C returns the index of the class with the highest membership probability: The parameters θ must be optimized to minimize the expected risk J (θ) defined by: where L is the per-example loss function quantifying the error between the predicted class C(X; θ) and the true class Y . However, as the true data-generating distribution p X,Y is unknown, the expected risk must be estimated by computing the expectation with respect to the empirical distribution p X,Y : Therefore, the training process consists in minimizing the empirical risk : 1 Ability to learn classes with complex boundaries (related to model complexity). This is known as the empirical risk minimization (ERM) principle [29]. Common choices for the function L are the hinge loss (defined for multiclass problems in [30]) used by support vector machines (SVMs), and the log loss or negative log-likelihood: widely used for classifiers based on artificial neural networks (ANNs) and for logistic regression. When L is the negative log-likelihood, the objective function J (θ) is the cross-entropy loss and the optimal set of parameters θ * minimizing J is the maximum likelihood estimator [31]. Usually, a regularization term is added to the empirical risk to penalize the model complexity in order to reduce overfitting. The boundaries between classes in the input space are called decision boundaries. Linear classifiers are classification algorithms for which the decision boundaries are defined by linear combinations of the features of X. Linear classifiers are appropriate when the classes are linearly separable in X , which means that the decision boundaries correspond to portions of hyperplanes. Linear classifiers include logistic regression [32,33,34], linear discriminant analysis (LDA [31]), and the linear support vector classifier (linear SVM [35,36]).
Many algorithms exist for nonlinear classification problems, each of them having its own advantages and drawbacks. As a kernel method, the linear SVM is extended to nonlinear classification problems using the kernel trick based on Mercer's theorem [37]. Artificial neural networks [38,39] (see [40] for a historical review) have become very popular due to their performances in numerous classification contests. Decision trees (e.g. CART algorithm [41]) and naive Bayes classifiers [42,43] are well-known for their interpretability. Other nonlinear classifiers include the k-nearest neighbors algorithm (kNN [44]), and quadratic discriminant analysis (QDA [31]). In [45], the most common classifiers are compared on eleven binary classification problems. Short reviews of classification algorithms can be found in [46,47].
Usually, combining several models to form a meta-estimator results in more robust predictions and reduces overfitting. This idea is used in ensemble methods such as bagging (or bootstrap aggregating) [48], feature bagging (or random subspace method) [49], stacking [31,50], boosting (including the well-known AdaBoost algorithm [51,52]), gradient boosting [53,54,55,56], and voting classifiers based on either a majority vote or a soft vote (technique known as ensemble averaging [57]). Random forests [58] combine bagging and feature bagging to build an ensemble of decision trees.

Classification for numerical simulations
Classification algorithms have recently found applications in numerical simulations, and more specifically for the selection of numerical models adapted to the context of the simulation. In this case, the class labels are used to identify the models.
Applications to turbulence modeling in computational fluid dynamics can be found in [7,10]. In large eddy simulations (LES, see [59]), the Navier-Stokes equations are filtered to avoid resolving small-scale turbulent structures whose effects are taken into account either by subgrid scale models (explicit LES closures) or via the dissipation induced by numerical schemes (implicit LES). In [7], sub-grid statistics obtained from direct numerical simulations enable training a fully-connected deep neural network to switch between different explicit LES closures at any point of the grid. This classifier is reused in [10], this time for switching between different numerical schemes in implicit LES. In both cases, the classifier is used to increase the accuracy of numerical predictions.
The idea of locally switching between different simulation strategies can also be found in [6] for the multiscale modeling of composite materials. In the multilevel finite-element method (FE 2 [60]), the quantities of interest at every integration point of the macroscopic finite-element mesh are given by a microscopic finite-element computation of an elementary cell representing the material's microstructure. The multi-fidelity surrogate model presented in [6] relies on two surrogate models replacing the microscopic finite-element model, namely a reduced-order model taken from [61] and an artificial neural network based regression model. At each integration point of the macroscopic mesh, the classifier (a fully-connected network) analyzes the effective strains and predicts whether the error of the regression model would be acceptable, enabling the selection of either the purely data-driven regression model or the more sophisticated physicsdriven ROM. This time, automatic model recommendation by a classifier is used to adapt the model complexity and reduce the computation time.
In [8,9], optimal classification trees (OCTs [62]) are used as model selectors in a data-driven physics-based digital twin of an unmanned aerial vehicle (UAV). The OCTs enable the update of the digital twin according to sensor data by selecting a model from a predefined model library. In this context, the training procedure for the classifier corresponds to an inverse problem. Indeed, training examples are generated by running simulations with all the models in the library and evaluating their predictions at the sensors' locations. Hence, for a given model y ∈ [[1; K]], the data x are obtained by means of numerical simulations performed with y. This corresponds to the forward mapping. The classifier must learn the inverse mapping giving y as a function of x. In this example, data labeling is straightforward: the label of a training example x is given by the index y of the model which was used to generate x. It is also noteworthy that generating training examples is not too expensive, because numerical simulations are performed with reduced-order models obtained by the Static-Condensation Reduced-Basis-Element method (SCRBE [63,64,65,66]). In this application, automatic model recommendation gives the UAV the ability to dynamically evaluate its flight capability and replan its mission accordingly.
Another example of classifier used to accelerate numerical simulations can be found in [4]. Contrary to [8,9], the data labeling procedure relies on the clustering of simulation data. In this framework, the model library is made of cluster-specific DEIM 2 [67] models that are faster than the high-fidelity model. The high-fidelity model computes a prediction u i for each input x i in the database {x i } 1≤i≤N , resulting in a dataset {u i } 1≤i≤N on which a clustering algorithm is applied. The predicted variable u is the discretization of a continuous field on a finite-element mesh, thus living in a high-dimensional space. To avoid the so-called curse of dimensionality [28], a DEIMbased feature selection technique is used before applying k-means clustering [68]. Alternatively, the clusters can be obtained with a variant of k-means using the DEIM residual as clustering criterion. Then, for a given training example x i , the class label y i is defined by the index of the cluster that u i is assigned to. In the exploitation phase, when dealing with test data, the best DEIM model is selected by a nearest neighbor classifier. The input data given to the classifier are either parameters of the problem or the variable u obtained at the previous time increment. A similar methodology is described in [5], where the concept of model library is termed model dictionary, which is the terminology adopted in this paper. The model dictionary is made of hyper-reduced-order models [69], and the input data {x i } 1≤i≤N are images of a mechanical experiment. The dimensionality of simulation data is reduced by Principal Component Analysis (PCA) before using k-means clustering. A convolutional neural network [70] is trained to return class labels without computing the intermediate variable u in order to avoid time-consuming operations. This classifier is an approximation of the true classifier K returning the correct label for any input x.

Definition of the classification problem
Notations: the j-th feature of a random vector X is the real-valued random variable denoted by X j . Its observations are denoted by x j , or x j i when indexing is necessary, for example when considering training data. When X is obtained by discretizing a random field on a mesh, the feature X j corresponds to the value taken by the random field at the j-th node. In the numerical application presented in this work, a random temperature field is considered. The spatial coordinates of the j-th node are stored in a vector ξ j ∈ R 3 . The categorical variable Y indicates which model should be used.
In this paper, input data {x i } 1≤i≤N correspond to several instances or variabilities of a physical field discretized on a mesh. Let N be the number of nodes in the mesh. If the physical field is scalar and defined at the nodes, then each observation x i is a vector of R N . For relatively small problems, N is in the order of 10 4 to 10 5 . For some industrial problems, N can be in the order of 10 6 to 10 8 . The dataset {x i } 1≤i≤N may come from experiments, numerical simulations, statistical models, or a combination of them, and contains from 10 2 to 10 4 observations. It is assumed that all features of all observations are known, contrary to some classification tasks in other disciplines encountering the problem of missing values. This assumption is clearly satisfied when data come from numerical simulations or statistical models. For experimental data, numerous techniques provide space-distributed measurements that can be projected onto the mesh, such as particle image velocimetry [71] in fluid dynamics, digital image correlation [72] and photoelastic experiments [73] in solid mechanics, and temperaturesensitive paints [74] measuring surface temperatures.
The framework considered in this paper is the same as in [11] for ROM-nets, where the input variabilities are supposed to be used for an uncertainty propagation study in a physics problem P, for which a high-fidelity model m HF is available. The physics problem P is a time-dependent problem. As the high-fidelity model is too computationally expensive, dictionary-based ROMnets have been introduced to reduce the computation time by means of a reduced-order model dictionary and a classifier playing the role of a model selector. The dictionary-based ROM-net is trained on the available dataset {x i } 1≤i≤N . For a given observation x i , the class label y i indicates the most appropriate model in the dictionary to be used for fast simulations with limited errors with respect to the high-fidelity model m HF . Class labels are obtained by the following data labeling procedure: • Step 1: for each observation x i in the dataset, use the high-fidelity model m HF to solve a simplified version P of the physics problem P (for example, the problem P can consist in solving P for a few time increments only). The primal solution of P computed for x i is denoted by u i . It consists of a collection {u n i } 1≤n≤nt of n t fields defined on the mesh, with n t being the number of time increments in problem P .
• Step 2: given {u i } 1≤i≤N , compute the dissimilarity matrix δ ∈ R N ×N with the following formula: with d Gr(∞,∞) being the Grassmann metric defined in [75]. The coefficient δ ij is a dissimilarity measure between x i and x j .
• Step 3: k-medoids clustering [76,77,78] is applied to the dissimilarity matrix δ. In this paper, we consider K = 4 clusters. The label is given by the index of the cluster containing u i .
. This dataset is split in a training set, a validation set and a test set with cardinalities 600, 200 and 200 respectively, enabling the supervised training and evaluation of a classifier C. For the sake of simplicity, the labeled data are renumbered so that the N train = 600 first input-output pairs {(x i , y i )} 1≤i≤N train form the training set on which the feature selection and data augmentation algorithms presented in this paper are trained.
In this work, the physics problem P is a temperature-dependent mechanical problem. The structure is made of an elasto-viscoplastic material whose behavior depends on the local value of the temperature field [79]. The random variable X is a random vector representing the evaluation of the random temperature field on a finite-element mesh containing N = 42445 nodes (see Figure 1). The structure is subjected to centrifugal forces and pressure loads. The random temperature fields are generated by a stochastic model described in [11], where ten fluctuation modes are randomly combined and superposed to a reference temperature field. The realizations of the random temperature field are continuous and always satisfy the heat equation. Modeling random fields as random combinations of deterministic spatial functions is quite common when studying stochastic partial differential equations [80,81,82], because a random field can be approximated by truncating its Karhunen-Loève expansion [83].
As already stated, the main contributions of this paper are a feature selection strategy and a data augmentation algorithm adapted to the specificities and difficulties of classification problems encountered when training dictionary-based ROM-nets. Concerning feature selection, the main focus is on the fast quantification of features redundancy by taking advantage of the type of input data. Concerning data augmentation, in addition to the constraints that have already been mentioned, it is likely that transforming an input example x i substantially modifies the intermediate variable u i , and thus the class label y i might no longer be relevant for the transformed input. Avoiding this situation is crucial to ensure that the augmented data are correctly labeled. Our algorithms are applicable under the assumptions that the random vector X derives from a random field whose realizations are continuous with probability one (sample path continuity, see Definition 2.1 in [84]) and belong to a convex domain X related to physics constraints. Lastly, a comparison of various classification algorithms is conducted to put into perspective the choice made in [11] to use an ensemble of deep neural networks trained with different architectures and loss functions. Remark 3.1. Another strategy would consist in using a regression algorithm for the classification task. Indeed, since our data labeling procedure is based on clustering, the classification problem could be replaced by a regression problem for the prediction of dissimilarities {δ(x,x k )} 1≤k≤K for x ∈ X , withx k being the medoid of the k-th cluster. Given these distances for a new observation x, the class label is obtained by taking the integer k ∈ [[1; K]] associated to the smallest dissimilarity δ(x,x k ). However, the data augmentation algorithm presented in this paper is not compatible with regression algorithms. For this reason, this paper focuses on classifiers rather than regressors.

.1 Feature selection based on mutual information
We recall that a projection π is a linear map satisfying π • π = π. It is entirely defined by its kernel and its image, which are complementary: given two complementary vector subspaces V 1 and V 2 , there is a unique projection π whose kernel is V 1 and whose image is V 2 , namely the projection onto V 2 along V 1 . For more details about projections, see [85], pages 385 to 388. Let us now give a formal definition of a feature selector : , the feature selector π S,B : V → V is the projection whose image is span ({e i } i∈S ) and whose kernel is span When the choice of the basis B is obvious, the notation π S,B is simply replaced by π S . In practice: Therefore, from a numerical point of view, one can interpret the feature selector as linear map π S : V → span ({e i } i∈S ), which enables reducing the size of the vector representing π S (x) for x ∈ V . In this way, applying a feature selector π S to a vector of R N consists in masking its features whose indexes are not in S, which gives a reduced vector in R |S| where |S| denotes the number of elements in S. Feature selection algorithms build the set S by searching for the most relevant features for the prediction of the output variable Y . For this purpose, the mutual information can be used to quantify the degree of the relationship between variables: Definition 4.2. (Mutual information [86], eq. 8.47, p. 251) Let Z 1 and Z 2 be two real-valued random variables with joint probability distribution p 1,2 and marginal distributions p 1 and p 2 . The mutual information I(Z 1 , Z 2 ) is defined by: The mutual information measures the mutual dependence between two random variables. Contrary to correlation coefficients, the information provided by this score function is not limited to linear dependence. The mutual information is nonnegative, and equals to zero if and only if the random variables are independent. Given Equation (2), replacing Z 1 by a feature X i of X and Z 2 by Y gives: The mutual information can be used to quantify the redundancy of a set of features S and its relevance for predicting Y : The minimum redundancy maximum relevance (mRMR) algorithm [21,22] builds the set S by maximizing D(S, Y ) − R(S), which is a combinatorial optimization problem. For this type of optimization problem, a brute-force search is intractable, because the number of solution candidates is too large. Instead, mRMR searches for a sub-optimal solution by following a greedy approach. First, the feature having the highest mutual information with the label variable Y is selected. Then, the algorithm follows an incremental procedure: given the set S m−1 obtained at iteration m − 1, form the set S m such that: This incremental procedure stops when m reaches the target number of features N f . A review of feature selection algorithms based on mutual information can be found in [87].

A geostatistical variant of mRMR feature selection
When training dictionary-based ROM-nets, the number of features of the random vector X scales with the number of nodes N in the mesh. In particular, the number of features is exactly N if X is the nodal representation of a scalar field. Hence, there are too many features to compute all redundancy terms I(X i , X j ). However, one can estimate the redundancy terms thanks to the proximities of the features on the mesh. Indeed, X is a regionalized variable: in our example, we recall that ξ i ∈ R 3 denotes the position of the i-th node in the mesh, and that the feature X i corresponds to the value taken by a random temperature field at ξ i . If two points ξ i and ξ j of the mesh are close to each other, the corresponding features X i and X j are likely to be correlated and thus redundant because of the smoothness of the temperature field. This idea is also valid when considering physical variables discretized in time.
In this paper, the random temperature field is modeled by a Gaussian random field [84] as in [11], which is a common and simple approach when modeling uncertainties on a physical field. As a consequence, X is a Gaussian random vector and the mutual information I(X i , X j ) has a simple formula involving the correlation coefficient: Property 4.5. (Mutual information of two correlated Gaussian random variables [86], eq. 8.56, p. 252) Let (X 1 , X 2 ) be a Gaussian random vector. The mutual information I(X 1 , X 2 ) reads: where ρ denotes the correlation between X 1 and X 2 .
This property implies that, for Gaussian random fields having isotropic correlation functions 3 ρ, the mutual information I(X i , X j ) only depends on the distance ||ξ i − ξ j || 2 . A wide variety of isotropic correlation functions are given in [84]. More generally, since Equation (14) is an increasing function of ρ 2 , any isotropic upper (resp. lower) bound of the squared correlation function gives an isotropic upper (resp. lower) bound of the mutual information.
For the example studied in this paper, Figure 2 shows that the mutual information I(X i , X j ) decreases as the corresponding distance ||ξ i − ξ j || 2 increases. Therefore, our feature selection algorithm builds a metamodelĨ replacing I(X i , X j ) by a function of the distance ||ξ i − ξ j || 2 , which drastically reduces the computational cost of mRMR algorithm for our particular problem. First of all, one must build a design of experiments (DOE) to select a few terms I(X i , X j ) to be computed exactly. The metamodelĨ is calibrated to fit the corresponding precomputed redundancy terms. Then, mRMR feature selection is applied by replacing I(X i , X j ) byĨ(||ξ i − ξ j || 2 ). The feature selection algorithm is described in Algorithm 1. We call this algorithm geostatistical mRMR, since geostatistics is the branch of statistics that deals with regonalized variables. A stopping criterion is added to the incremental procedure used in mRMR, enabling an automatic selection of the number of features to be kept for the classification task: the algorithm stops when the value of arg max A condition on the mutual information I(X i , Y ) can also be added to avoid selecting quasi-irrelevant features. For stage 1 of Algorithm 1, computing all the terms ξ j 1 − ξ j 2 2 of the matrix of pairwise mesh nodes distances is not necessary: only a few lines of this matrix corresponding to randomly selected nodes are evaluated, which is sufficient to build the DOE. In other words, one computes the distances between a few nodes and all the mesh nodes.
Remark 4.6. A parallel can be drawn between our feature selection strategy and hyper-reduction methods [69,88,89,90] used to accelerate complex nonlinear problems in physics (see [91] for design optimization and [92] for large-scale simulations). Hyper-reduction methods aim at finding a reduced set of integration points in the finite-element mesh that is sufficient to predict the behavior of the physical system. The constitutive equations are solved on this reduced integration domain only, while the values of quantities of interest at the remaining integration points can be recovered with the Gappy-POD [93]. In short, hyper-reduced solvers make predictions from a reduced number of points in a mesh, like the classifiers used in this paper do when combined with the geostatistical mRMR. Although the objectives are different, both hyper-reduction and geostatistical mRMR feature selection benefit from the properties of physics data to reduce the complexity of numerical tasks.

Numerical results
The red curve on Figure 2 corresponds to the metamodel estimating redundancy terms. In this example, we choose: Select distance values r j .

5:
Compute the mutual information I(X i , X j ) for each pair selected in Stage 1.

6:
Train a metamodelĨ such that I(X i , X j ) ≈Ĩ ξ i − ξ j 2 .  where H is the Heaviside step function and I ∞ , γ 1 , γ 2 , r 1 , r 2 , α 1 , α 2 are calibration parameters that are adjusted manually. In the DOE, the step between distances r j is smaller for small distances, in order to better capture the evolution of the mutual information in its high gradient regime. The number n j of pairs of nodes separated by a distance of r j selected in the DOE also depends on r j : as higher variances were expected for small distances, n j decreases when r j increases. In total, 749 terms I(X i , X j ) are computed, which takes 5.12 seconds using Scikitlearn [94]. Building the DOE takes only 0.33 seconds. Then, the greedy procedure takes 303 seconds and selects 87 features among the 42445 original ones. The first iteration is the longest one with 276 seconds, because it includes the computation of all the relevance terms I(X i , Y ). As a comparison, the original mRMR algorithm takes 6469 seconds to compute 7 iterations only. We did not let mRMR algorithm go further, since the per-iteration computation time grows with the iteration number. For a fair comparison, our implementations of mRMR and stages 3 and 4 of the geostatistical mRMR are the same except that redundancy terms are evaluated with Scikit-learn for mRMR and with the functionĨ for the geostatistical mRMR. Table 1 compares the relevance D(S, Y ), the true redundancy R(S), the approximate redun-dancyR(S) estimated withĨ, the true cost function D(S, Y ) − R(S) and the approximate cost function D(S, Y ) −R(S) for three different feature selection strategies: • the geostatistical mRMR feature selection (Algorithm 1), selecting a set S * of features; • a univariate filter algorithm selecting the features with the highest mutual information (MI) scores I(X i , Y ). This algorithms finds a set S M I maximizing the relevance for a given cardinality; • a purely geometric feature selection algorithm, randomly selecting the first feature and adding features in a greedy manner so that the distance to the closest point ξ i , i ∈ S m is maximized. This algorithm tends to select a set S G of well-distributed features in order to get a low redundancy for a given cardinality.
Since the geostatistical mRMR automatically selected 87 features, the two other approaches are applied with |S G | = |S M I | = 87 as a target. Table 1 shows that the relevance of the set S * selected by our algorithm is in the same order of magnitude as the relevance of the set S M I . Its redundancy is in the same order of magnitude as the redundancy of the set S G . These results show that the geostatistical mRMR algorithm does have the desired behavior: it selects a subset of features S * with high relevance and low redundancy. Figure 3 shows the features selected by the three different algorithms. The classification accuracies of several classifiers using the reduced features S * are given in the last section of the article.  Remark 4.7. The geometric feature selection algorithm gives rather good results in terms of the cost function, but it does not mean that it is an appropriate approach. Indeed, one can see that the relevance of S G is very low, since this algorithm does not use any information concerning the classification problem.  [95], p. 10) Let V be a real vector space. A non-empty set S ⊂ V is convex if: Definition 5.2. (Convex combination [95], p. 11) Let {x i } 1≤i≤n be a finite set of elements of a real vector space V . A convex combination of {x i } 1≤i≤n is a vector x ∈ V such that: Proof. Let z ∈ E (L(S)). Following the definition of a convex hull, there exists n ∈ N * such that: For all i ∈ [[1; n]], as w i ∈ L(S), there exists v i ∈ S such that w i = L(v i ). By linearity of L: so E (L(S)) ⊂ L (E(S)). The other inclusion can be shown using exactly the same arguments. Thus: L (E(S)) = E (L(S)).
This property has a very simple yet important consequence for the data augmentation algorithm presented in this paper: Property 5.5. Let V and W be two real vector spaces, and let L : V → W be a linear map. Let S be a non-empty set included in V . Then, for all x ∈ V : Proof. By contraposition, x ∈ E(S) ⇒ L(x) ∈ L (E(S)) = E (L(S)).
Our data augmentation strategy uses this property in the particular case where the linear map is a projection. As a reminder, the notation K stands for the true classifier assigning any input For any subset S k of C k with cardinality |S k |,Â S k ∈ R dim(V )×|S k | denotes the matrix whose columns contain the coordinates of the elements of S k . The matrix denoted by A S k is obtained by adding a row of ones below the matrixÂ S k , giving a matrix of size (1 + dim(V )) × |S k |.
Theorem 5.7. (Pure set characterization) Let S k be a subset of C k with cardinality |S k |. The set S k is pure in S if and only if for all x in S \ C k the linear system: has no nonnegative solution w ∈ R |S k | + . Proof. Let x ∈ S \ C k . Equation (23) has no nonnegative solution if and only if: which ends the proof.
Corollary 5.8. (Pure set testing) Let S k be a subset of C k with cardinality |S k |, and let L : V → W be a linear map, where W is a finite-dimensional real vector space. If for all x in S \ C k the linear system: (26) has no nonnegative solution in R |S k | + , then S k is pure in S. Proof. Equation (26) characterizes the purity of L(S k ) in L(S) (Theorem 5.7), which implies that S k is pure in S (Property 5.5). Figure 4 illustrates the concept of pure sets. On this figure, the set C 1 is made of all the elements represented by dots, while the crosses form the set C 2 = S \ C 1 . On the left, the subset formed by the six black dots is pure since its convex hull delimited by dashed lines contains only dots. The subset made of the six black dots on the right-hand side of the figure is not pure because of the presence of a cross in its convex hull. Equation (23) has a nonnegative solution when using the coordinates of this cross in its right-hand side.

The data augmentation algorithm
The objective is to generate new data points x ∈ X in a given class y ∈ [[1; K]] from the preexisting observations in that class. To this end, one must apply class-preserving transformations on the training examples. New examples can be created by taking convex combinations of some subsets of the training set, for example. One way of controlling the risk that newly generated examples have wrong labels is to take convex combinations of subsets only if they are pure. Indeed, if the k-th class C k contains a set S k that is pure in the training set, one can expect that the probability P(Y = k | X ∈ E (S k )) is high enough to get new examples of class C k by drawing samples in E (S k ). In addition, the third Hadamard well-posedness condition states that the solution of a physics equation changes continuously with respect to the parameters of the problem. In the neighborhood of a point x 0 belonging to a pure set S k , the primal solution u stays in the neighborhood of the solution u 0 obtained with x 0 and is thus likely to have the same label. Hence, the objective of our algorithm is to find pure sets in the training set in order to generate new examples by convex combinations with a limited risk of getting incorrectly labeled examples. The pure sets detected by the algorithm are listed in a matrix S such that S[k, i] contains the indices of the training points forming the i-th pure set of the k-th class. The pure sets are grown from different starting points or seeds in the training set by iteratively adding the seeds' nearest neighbors in terms of the precomputed dissimilarity measure δ used for clustering in the data labeling procedure. The growth stops before losing the purity of the subsets.
However, checking the purity in the high-dimensional input space can cause difficulties, even when training the data augmentation algorithm after a first dimensionality reduction like in this paper. For this reason, the algorithm checks the purity after having applied a feature selector π S with a small random subset of features S containing d features. Let us apply Property 5.5 with V = W being the input vector space containing X and with the linear map L being the feature selector π S . As Property 5.5 states, if no point of π S ({x m } 1≤m≤N train \ C k ) belongs to the convex hull of π S {x m } m∈S[k,i] , then the set E {x m } m∈S[k,i] does not contain any point labeled with k = k. Since a set can lose its purity after projection, the algorithms tries p max random feature selectors π S before considering that the set is not pure. In practice, the purity of π S {x m } m∈S [k,i] in π S ({x m } 1≤m≤N train ) is numerically tested by solving a nonnegative least squares (NNLS [96]) problem. If for all points x ∈ {x m } 1≤m≤N train \ C k the inequality: is satisfied withπ S (x) = (π S (x) T 1) T and with ε DA being the tolerance of the data augmentation algorithm, then Corollary 5.8 implies that {x m } m∈S[k,i] is pure in {x m } 1≤m≤N train . Algorithm 3 describes the data augmentation algorithm. It calls Algorithm 2 to find n well-distribued seeds per class before growing pure sets. It is noteworthy that using few pure sets to generate many examples would increase the distribution gap [97] between augmented data and original data. To avoid this issue, one had better use many well-distributed seeds to distribute data augmentation efforts between the pure sets and thus limit the divergence between the augmented distribution and the true data-generating distribution.
Remark 5.9. Realizations of the random variable X belong to a convex domain X related to physics constraints. When considering surface random temperature fields defined on the boundaries of a solid, X is a hypercube consisting of all the fields taking values between zero Kelvin degree and the material's melting point. These random fields can be used as Dirichlet boundary conditions for the nonlinear heat equation. The assumption of a linear thermal behavior is added when considering three-dimensional random temperature fields defined inside the solid, so that the set X remains convex when adding the constraint that the random field must satisfy the heat equation. More generally, convex combinations respect physics constraints defined by linear operators, such as linear partial differential equations and Dirichlet, Neumann and Robin boundary conditions.

Numerical results
Linear discriminant analysis (LDA), commonly used for classification tasks, can also be used for supervised dimensionality reduction by projecting the data onto the subspace maximizing the between-class variance, as explained in [31]. For the classification problem presented in this paper, the training data are visualized in the two-dimensional subspace obtained by LDA in Figure 5. Although this subspace is the one that best separates the classes, one can see that the training examples do not form well-separated groups. For this reason, testing the purity of subsets of training data before generating new examples by convex combinations is necessary to reduce the risk of getting incorrectly labeled augmented data. The data augmentation algorithm finds about 60 pure sets per class with an average population of 5 training examples, using random subspaces of dimension 5 to test the purity. Note that two pure sets are merged only when one is included in the other, since the union of two pure sets is not always pure. The computation time for the data augmentation training phase is 40 minutes. Once pure sets have been found, one can generate as many augmented examples as necessary. Generating 5400 examples to multiply the size of the training set by 10 takes less than a second. Among the augmented data, 400 examples are taken for the evaluation of the Algorithm 2 Seeds selection for data augmentation. Note: all the dissimilarities have already been computed in the data labeling procedure.
Input: training set {(x i , y i )} 1≤i≤N train , class label k, class centerx k , dissimilarity matrix δ, target number of seeds n, preselection parameters (ε 1 , ε 2 ) ∈ [0; 1] 2 . Output: List l k of n indices of seeds candidates for the k-th class. Append j to l k . 11: end for 12: return l k Figure 5: Data visualization in the 2D subspace maximizing the separation between classes (supervised linear dimensionality reduction using LDA). data augmentation algorithm. The data labeling procedure involving numerical simulations is applied for these 400 examples in order to estimate the percentage of incorrectly labeled data. It turns out that none of these examples is incorrectly labeled, which validates the algorithm for our problem. The benefits of data augmentation for the classification task are evaluated in the final section.
6 Validation of our feature selection and data augmentation algorithms 6 Apply Algorithm 2 to get the list l k of n indices of seeds candidates. purposes, each classifier is tested twice: once in combination with the geostatistical mRMR and once with principal component analysis (PCA) with 10 modes. Since the random temperature fields derive from a Gaussian random field involving only 10 modes, the database obtained after applying PCA contains all the information of the original data. Each combination of one of the 14 classifiers with PCA or feature selection is trained twice: once on the true training set containing N train = 600 examples, and once on the augmented training set made of 6000 examples.
All the classifiers are trained with Scikit-learn [94], except multilayer perceptrons (MLPs) and radial basis function networks (RBFNs) which are trained with PyTorch [98]. We train the RBFNs in a fully supervised manner, which means that the parameters of the radial basis functions are learnt by gradient descent like the weights of the fully-connected layers. In addition, we use only one hidden layer for RBFNs, since these artificial neural networks generally have shallow architectures, as explained in [99]. Deeper architectures have been tested for MLPs. Scikit-learn's MLP classifier has also been tested; it is called simple MLP in this paper, because its architecture is only made of fully-connected layers and does not include dropout [100] nor batch normalization [101]. All the classifiers based on artificial neural networks are trained with Tikhonov regularization and early stopping [19]. Logistic regression is trained with elastic net regularization [102]. Kernels used for support vector machines (SVMs) are obtained by combining several polynomial kernels with different hyperparameters. Kernel design could be optimized using multiple kernel learning algorithms [103], but we simply build our kernels by evaluating different combinations on the validation set, just as when we look for a good architecture for artificial neural networks.
The classification accuracies on test data are given in Table 2. Of course, this ranking is specific to the classification problem presented in this paper, no general conclusion can be drawn from this particular numerical application. On this classification problem, when using augmented data in the training phase, the highest test accuracy reached with linear classifiers is 43.5%, obtained with the linear SVM combined with PCA. The fact that k-nearest neighbors classifiers barely exceed 50.0% of accuracy on this problem is related to an observation that was made in [11]: there is no simple correlation between the Euclidean distance and the physicsinformed dissimilarity measure used in dictionary-based ROM-nets. MLPs get the best results, reaching 87.0% of accuracy when combined with our data augmentation and feature selection algorithms. Interestingly, quadratic discriminant analysis (QDA) gives excellent results while having no hyperparameter to tune, contrary to the two other families of classifiers obtaining the best results, namely MLPs and multiple kernel SVMs. This makes QDA the best compromise between accuracy and training complexity for this specific classification task.
Although PCA perfectly describes the input data in this example, the geostatistical mRMR feature selection algorithm enables reaching higher accuracies with some classifiers. Not only it behaves as the original mRMR when selecting features, but it also gives satisfying results when combined with a classifier. Concerning data augmentation, Table 2 shows that our algorithm significantly improves classification results. The accuracy gain due to data augmentation is 4.98% on average and ranges from −2.5% to 10.5%, increasing the accuracy in 25 cases out of 28.

How to further improve classification performances?
Ensemble methods can be used to reduce overfitting and increase the accuracy on test data. In addition, it enables recycling different variants of a classifier that the user has trained for different hyperparameters. Using ensemble averaging with classifiers trained on the augmented dataset with feature selection, we manage to combine 6 MLPs with different architectures to reach an accuracy of 89.0%. When stacking these MLPs with a ridge logistic regression analyzing the predicted membership probabilities, we get an accuracy of 90.0%. In addition to ensemble learning methods, one can also use random noise injection to increase noise robustness, as

Conclusion
Classification algorithms are used in computational physics for automatic model recommendation. Such modeling strategies enable the reduction of the computation time, or the selection between models with different physics when one wants to improve the accuracy of numerical predictions. This article deals with the specificities of the classification problems encountered in computational physics, and more particularly for dictionary-based ROM-nets. These classification problems generally have the three following issues: the lack of training data, their high dimensionality, and the non-applicability of common data augmentation techniques to physics data. To tackle these difficulties, two algorithms are proposed. The first one is a geostatistical variant of the mRMR feature selection algorithm, enabling the identification of a reduced set of relevant but non-redundant features for high-dimensional regionalized variables. The second one is a data augmentation algorithm controlling the risk of generating new examples with wrong labels by finding pure subsets in the training set. The performances and benefits of these algorithms are illustrated on a classification problem for which 14 classifiers are evaluated.