We assume that a sufficiently large database is available, where a physical property of interest and a number of associated ruling primitive variables or observables are stored.

We introduce and test two machine learning approaches to discover possible groups or combinations of primitive variables: The first approach is based on regression models whereas the second on classification models.

The variable group (here referred to as the new effective good variable) can be considered as successfully found, when the physical property of interest is characterized by the following effective invariant behaviour: In the first method, invariance of the group implies invariance of the property up to a given accuracy; in the other method, upon partition of the physical property values into two or more classes, invariance of the group implies invariance of the class.

For the sake of illustration, the two methods are successfully applied to two popular empirical correlations describing the convective heat transfer phenomenon and to the Newton's law of universal gravitation.

Keywords Machine Learning in Physics, Primitive Variable Analysis, Physical Property Invariance, Feature grouping

§ INTRODUCTION

Theoretical modelling and numerical simulations have become invaluable tools for the current scientific and technological advancement across various fields [1].

In their essence, models make use of a mathematical description of the physical laws by establishing relationships between physical variables.

The latter variables have the key role of describing in a complete and non-redundant manner a system of interest.

As such, their correct identification is never a trivial task, especially in systems with little prior knowledge [2, 3, 4].

Often, in order to effectively model complex systems, it is favourable to search for a more convenient description with a reduced number of effective variables [5, 6].

In other words, the system behaviour is not described by the directly accessible physical variables, but rather by some groups or combinations.

In this respect, back in the late 19th century, the Buckingham theorem was introduced [7, 8, 9].

In the latter approach, based on dimensional analysis, it is possible to mix the primitive physical variables, thus creating fewer dimensionless numbers which are effectively relevant.

Over the years, sophisticated methodologies have been developed also based on machine learning, data mining and other data-driven approaches [10, 11, 12].

Some other approaches aim to unlock the discovery of symbolic expressions accurately matching data derived from an unknown function.

This problem has been tackled with a number of methods [13, 14], including sparse regression [15, 16, 17, 18], genetic algorithms [19, 20, 21] and physics-inspired algorithms like AI Feynman, introduced by Udrescu & Tegmark [22].

Inspired by the latter, some authors of this work presented a multi-objective optimization procedure for reducing the set of composition-based material descriptors by optimally mixing them in power combination form.

This resulted in improved classification performances, as demonstrated in a case study focused on superconductors [23].

Moreover, the same procedure was applied by Bonke et al. [24] to identify an effective reduced set of variables in a micelle-based photocatalytic system for solar fuels production.

Specifically, the algorithm allowed to analytically mix five physical primitive variables into two synthetic features for the optimal binary separation of the experimental samples according to their performances.

The practical benefit of such approach lies in its capability to replace an experimental sample achieving high performance  with an alternative combination of its elemental components.

This gives the possibility to reduce the most expensive ones re-balancing the others, at no cost on the overall performance.

In this work we present a general and automated methodology, suitable for regression tasks, able to identify groups and/or group sets of variables in the power form $x_i^{\alpha_1} x_j^{\alpha_2} \cdots x_m^{\alpha_p}$.

We show its effectiveness on popular thermo-fluid dynamics correlations (i.e. Dittus-Boelter and Gnielinski).

Furthermore, we also demonstrate that the approach can be easily extended to more general functional form, proving its ability on the Newton's law of universal gravitation.

Also, we employ the optimal feature mixing procedure in ref. [23] for classification tasks, showing its successful application on the former two problems of this study.

The paper is organized in sections as follows.

Section <ref> illustrates the procedure for generating the three datasets employed in this study, based on the three functional forms analyzed here (i.e., Dittus-Boelter correlation, Gnielinski correlation, Newton's law of universal gravitation).

Furthermore, it presents the methodology for searching good variables in regression models and describes the implementation of classification models for optimal variable mixing towards class separation.

This section is further subdivided into specific techniques, including the consideration of single and multiple invariant groups in power form, as well as a broader exploration into non-power forms.

Section 3 presents numerical examples and discussions related to our study.

We illustrate the application of our methodology on selected cases, including the analysis of the popular Dittus-Boelter and  Gnielinski correlations (each examined from both regression and classification perspectives in power group forms), and Newton's law of universal gravitation (only regression on generic group forms).

Here, we also provide a comprehensive evaluation of our approach.

Finally, in Section 4, we draw conclusions providing an overview of the main contributions.

Overview of the protocol used to detect possible symmetries of a target property of interest with respect to its input variables, utilizing only data and ignoring the analytical functional dependence. Two distinct methodologies are presented: the former for identifying, in regression tasks, invariant groups in the form $x_i^{\alpha_1} x_j^{\alpha_2} \cdots x_m^{\alpha_p}$, among others; the latter for identifying, in classification tasks, one or several mixed features as power combinations of the input variables to achieve an optimal class separation.

§ METHODS

§.§ Datasets creation

We first create three datasets from physical models to be used to train an test Machine Learning (ML) tools, aiming to predict physical properties of interest.

More specifically, the first two datasets are generated using popular thermo-fluid dynamic correlations (i.e. Dittus-Boelter and Gnielinski correlations) with the target quantity being the Nusselt number.

Those datasets are created starting from the physical properties of 16 real liquids at ambient temperature and pressure [25], and evaluating the kinematic viscosity $\nu$ and the thermal diffusivity $\kappa$ of each liquid from tabulated data (see Supplementary Note 1).

Each fluid accounts for 500 value combinations of the flow speed $u$, the hydraulic diameter $D$, and the friction factor $f$ (the latter only requested for the Gnielinski correlation).

The values for these three variables are randomly chosen in the ranges $[0.1,1]$, $[0.01,0.1]$, and $[0.02,0.09]$, respectively.

This leads to a total of 8,000 samples for each dataset.

The Nusselt number is thus calculated according to the corresponding equations.

Finally, emulating what is typically experienced in experimental measurements, noise is added on top of the correct target values, sampling 8,000 points $\eta_{i}$ from a Gaussian distribution with mean $\mu=0$ and standard deviation $\sigma=0.1$; the target value of each entry $\mathrm{Nu}_{i}$ is thus multiplied by $(1+\eta_{i})$.

The third dataset is created on the basis of the Newton's law of universal gravitation.

It is constructed in a similar fashion, by generating random values within the range $[10^{16},~10^{18}]$ for the masses $m_{1}$ and $m_{2}$, and within the ranges $[-10^{12},-10^{10}]$ and $[10^{10},10^{12}]$ for the spatial coordinates $x_{1}, y_{1}, z_{1}, x_{2}, y_{2}, z_{2}$.

Here, the target value is the gravitational force $F_\mathrm{g}$ which is computed with some noise added following the same procedure as above.

§.§ Searching for good variables by regression models

§.§.§ Single invariant group in power form

Let $f(\mathbf{x}) = f(x_{1},\dots,x_{n})$ denote a function depending on $n$ physical variables $x_{1},\dots, x_{n}$.

If $f(\mathbf{x})$ is invariant with respect to a group of $p$
variables in the form $(x_{i}^{\alpha_1}x_{j}^{\alpha_2}\cdots x_{m}^{\alpha_p})$, $\alpha_1, \alpha_2, \dots, \alpha_p \in\mathbb{R}$, when this group is kept at a constant value $\overline{c}$, regardless of the value of each component $x_i, x_j, \dots, x_m$, $f(\mathbf{x})$ does not change.

Stated differently, the above invariance condition requires:

\begin{equation}
\label{eq2}
\alpha_1\ln(x_i) + \alpha_2\ln(x_j) + \cdots + \alpha_p\ln(x_m) = c
\end{equation}

($c= \ln(\overline{c})$), which upon differentiation yields:

\begin{equation}
\label{eq3}
\alpha_1\frac{\mathrm{d}x_{i}}{x_{i}} + \alpha_2\frac{\mathrm{d}x_{j}}{x_{j}} + \cdots + \alpha_p\frac{\mathrm{d}x_{m}}{x_{m}} = 0.
\end{equation}

We can thus construct the $(1 \times p)$ matrix $\mathbf{B}$ at a generic point $\overline{\mathbf{x}_{0}} = (x_{i,0}, x_{j,0}, \dots, x_{m,0})$, as

\begin{equation}
\label{eq4}
\mathbf{B} = \left[ \frac{\alpha_1}{x_{i,0}}, \frac{\alpha_2}{x_{j,0}}, \cdots, \frac{\alpha_p}{x_{m,0}} \right],
\end{equation}

Let $\mathbf{K}$ be a matrix whose columns form an ortho-normal basis in the null space (or kernel) of $\mathbf{B}$.

Hence, the condition of invariance of $f(\mathbf{x})$ with respect to the group $(x_{i}^{\alpha_1}x_{j}^{\alpha_2}\cdots x_{m}^{\alpha_p})$ can be recast into an orthogonality condition in the configuration space between the normalized gradient of $f(\mathbf{x})$ and each column of K, namely

\begin{equation}
\label{eq5}
\left({\nabla \Tilde{f}} \right)_{\overline{\mathbf{x}_{0}}}   \cdot  {\mathbf{K}} = 0
\end{equation}

\begin{equation}
\label{eq6}
{\nabla \Tilde{f}} =\frac{\nabla f}{\mathrm{norm}(\nabla f)}
%     \hspace{.5cm}
%    \Tilde{\mathbf{K}}^\top_i = \frac{{\mathbf{K}}^\top_i}{\mathrm{norm}({\mathbf{K}}^\top_i)}.
\end{equation}

Coefficients $\alpha_1, \alpha_2, \dots, \alpha_p$ can be conveniently normalized thus yielding the following algebraic system:

\begin{equation}
\label{eq7}
\begin{cases}
%\left( \frac{\partial \Tilde{f}}{\partial x_{i}} %\right)_{\overline{\mathbf{x}_{0}}} \Tilde{K}_{11} + \left( %\frac{\partial \Tilde{f}}{\partial %x_{j}}\right)_{\overline{\mathbf{x}_{0}}} \Tilde{K}_{21} + %\cdots + \left( \frac{\partial \Tilde{f}}{\partial %x_{m}}\right)_{\overline{\mathbf{x}_{0}}} \Tilde{K}_{\Phi1} = 0 %\\
%\hspace{.5cm}\vdots \\
%\left( \frac{\partial \Tilde{f}}{\partial x_{i}} %\right)_{\overline{\mathbf{x}_{0}}} \Tilde{K}_{1(\Phi-r)} + %\left( \frac{\partial \Tilde{f}}{\partial %x_{j}}\right)_{\overline{\mathbf{x}_{0}}} \Tilde{K}_{2(\Phi-r)} %+ \cdots + \left( \frac{\partial \Tilde{f}}{\partial %x_{m}}\right)_{\overline{\mathbf{x}_{0}}} \Tilde{K}_{\Phi(\Phi-%r)} = 0 \\
\left({\nabla \Tilde{f}} \right)_{\overline{\mathbf{x}_{0}}}   \cdot  {\mathbf{K}} = 0
\\
\alpha_1^{2}+\alpha_2^{2}+ \cdots +\alpha_p^{2}=1
\end{cases}
\end{equation}

If the non-linear system in Eq. (<ref>) is satisfied for the same exponents ($\overline{\alpha_1}$, $\overline{\alpha_2}$, …, $\overline{\alpha_p}$) over all the domains of the features ($x_{i}$, $x_{j}$, …, $x_{m}$), the group $x_{i}^{\overline{\alpha_1}}x_{j}^{\overline{\alpha_2}}\cdots x_{m}^{\overline{\alpha_p}}$ does represent an intrinsic variable and $f(\mathbf{x})$ is invariant with respect to that group.

§.§.§ Multiple concurrent invariant groups in power form

Clearly, the above procedure can be extended to a function $f(\mathbf{x})$ being invariant with respect to a higher number of feature groups, with each group even sharing some of the primitive physical variables.

Without loosing generality, and for the sake of simplicity, we limit this generalized description to sets consisting of two concurrent groups.

Let $f(\mathbf{x}) = f(x_{1},\dots,x_{n}) = f(x_{1},\dots,x_{i}^{\alpha_1}x_{j}^{\alpha_2}\cdots x_{m}^{\alpha_p}, x_{j}^{\beta_1}x_{k}^{\beta_2}\cdots x_{r}^{\beta_q}, \dots, x_{n})$
denote a function where two groups share a generic primitive variable $x_j$.

In this case, the above procedure applied to individual groups $x_{i}^{\alpha_1}x_{j}^{\alpha_2} \cdots x_{m}^{\alpha_p}$ and $x_{j}^{\beta_1}x_{k}^{\beta_2} \cdots x_{r}^{\beta_q}$ independently is not suitable any longer and requires a generalization as discussed below.

To investigate the invariance with respect to both groups, we now construct the $(2\times l)$
matrix $\mathbf{B}$:

\begin{equation}
\label{eq8}
\mathbf{B} =
\begin{bmatrix}
\frac{\alpha_1}{x_{i,0}} & \frac{\alpha_2}{x_{j,0}} & 0 & \cdots & \frac{\alpha_p}{x_{m,0}} & 0 \\
0 & \frac{\beta_1}{x_{j,0}} & \frac{\beta_2}{x_{k,0}} & \cdots & 0 & \frac{\beta_q}{x_{r,0}}
\end{bmatrix},
\end{equation}

with $ \max(p, q) \le l \le p+q$, $\mathbf{K}$ being a matrix whose columns represent an ortho-normal basis of the null space of $\mathbf{B}$.

As done above, the condition of invariance requires that the normalized gradient of $f(\mathbf{x})$ is orthogonal to each column of the kernel matrix K (see also Eq. (<ref>)).

Adding the normalization condition of the coefficients $\alpha_1, \alpha_2, \dots, \alpha_p$ and $\beta_1, \beta_2, \dots, \beta_q$, we obtain:

\begin{equation}
\label{eq10}
\begin{cases}
%        \frac{\partial \Tilde{f}}{\partial x_{i}} \Tilde{K}_{11} + \frac{\partial \Tilde{f}}{\partial x_{j}} \Tilde{K}_{21} + \frac{\partial \Tilde{f}}{\partial x_{k}} \Tilde{K}_{31} + \dots + \frac{\partial \Tilde{f}}{\partial x_{m}} \Tilde{K}_{(\Phi-1)1} + \frac{\partial \Tilde{f}}{\partial x_{p}} \Tilde{K}_{\Phi1} = 0 \\
%        \\
%        \hspace{.5cm}\vdots \\
%        \\
%        \frac{\partial \Tilde{f}}{\partial x_{i}} \Tilde{K}_{1(\Phi-r)} + \frac{\partial \Tilde{f}}{\partial x_{j}} \Tilde{K}_{2(\Phi-r)} + \frac{\partial \Tilde{f}}{\partial x_{k}} \Tilde{K}_{3(\Phi-r)} + \dots + \frac{\partial \Tilde{f}}{\partial x_{m}} \Tilde{K}_{(\Phi-1)(\Phi-r)} + \frac{\partial \Tilde{f}}{\partial x_{p}} \Tilde{K}_{\Phi(\Phi-r)} = 0 \\
\left({\nabla \Tilde{f}} \right)_{\overline{\mathbf{x}_{0}}}   \cdot  {\mathbf{K}} = 0
\\
\alpha_1^{2}+\alpha_2^{2}+\dots+\alpha_p^{2}=1
\\
\beta_1^{2}+\beta_2^{2}+\dots+\beta_q^{2}=1
\end{cases}
\end{equation}

with the partial derivatives being evaluated in $\overline{\mathbf{x}_{0}}$.

If the non-linear system in Eq. (<ref>) is satisfied for the same exponents ($\overline{\alpha_1}$, …, $\overline{\alpha_p}$, $\overline{\beta_1} \dots, \overline{\beta_q}$) in the entire feature domain ($x_{i}$, $x_{j}$, …, $x_{r}$), the two groups $x_{i}^{\overline{\alpha_1}}x_{j}^{\overline{\alpha_2}} \cdots x_{m}^{\overline{\alpha_p}}$ and $x_{j}^{\overline{\beta_1}}x_{k}^{\overline{\beta_2}} \cdots x_{r}^{\overline{\beta_q}}$ are intrinsic variables as the function $f(\mathbf{x})$ is invariant with respect to both groups.

It is worth stressing that more general cases with three or more concurrent invariant groups will imply additional rows for the matrix $\mathbf{B}$.

Furthermore, in general, the non-linear system (<ref>) is not necessarily closed (see subsection <ref>).

§.§.§ Further generalization to non power forms

The invariant group/set identification procedure, introduced for groups in the power form, can be generalized to other functional relationships.

For the sake of illustration, in the following we restrict to a single invariant group.

Let $f(\mathbf{x}) = f(x_{1},\dots,x_{n})$ denote a function invariant with respect to a group involving $p$ variables according to a generic functional dependence $g(x_i, x_j, \dots, x_m)$.

When such group is a constant $\overline{c}$, even varying its components $x_i, x_j, \dots, x_m$ separately, $f(\mathbf{x})$ does not change.

This yields:

\begin{equation}
\label{eq11}
g(\underbrace{x_{i}, x_{j}, \dots, x_{m}}_{p}) = \overline{c}
\end{equation}

which translates into:

\begin{equation}
\label{eq12}
\frac{\partial g}{\partial x_{i}} \mathrm{d}x_{i} + \frac{\partial g}{\partial x_{j}} \mathrm{d}x_{j} + \dots + \frac{\partial g}{\partial x_{m}} \mathrm{d}x_{m} = 0.
\end{equation}

We can thus construct the $(1\times p)$ matrix $\mathbf{B}$ at a generic point $\overline{\mathbf{x}_{0}} = (x_{i,0}, x_{j,0}, \dots, x_{m,0})$, as

\begin{equation}\label{eq13}
\mathbf{B} = \left[ \frac{\partial g}{\partial x_{i}}, \frac{\partial g}{\partial x_{j}}, \dots, \frac{\partial g}{\partial x_{m}} \right]
\end{equation}

Let $\mathbf{K}$ be a matrix whose colomns represent an orthonormal basis of the null space of $\mathbf{B}$.

Applying the invariance condition of Eq. (<ref>) leads to a non-linear system analogous to Eq. (<ref>).

Notably, the procedure illustrated here can be further extended to sets of groups in any functional form, following a reasoning similar to that applied in the above subsection <ref>.

§.§.§ Regression model and procedure implementation

In this study, we assume that a database of physical data is available.

Each sample in our database consist of $n$ features $(x_1, \dots, x_n)$ corresponding to a target quantity $f(\mathbf{x})$.

We further assume that the target quantity possibly depends on some invariant groups expressed for instance in the power form.

We thus try to detect such invariance following the above procedure.

As already described, this requires the evaluation of the gradient $\nabla f(\mathbf{x})$.

To this end, $f(\mathbf{x})$ can be conveniently approximated with a Deep Neural Network (DNN), allowing the computation of that gradient by means of automatic differentiation.

All the DNNs of this study are trained and validated over the 85% of the corresponding databases – of which the 85% is used for the training and the remaining 15% for the validation – and tested over the remaining 15%.

The DNN structure, as reported in Supplementary Note 2, is used for all the examples of this study.

Upon network training and validation, automatic differentiation is utilized to calculate the gradient of the function at a specific point $\mathbf{x}_{0}$ within its domain.

Here, the variables in the examined group are randomly chosen within their domains, while the values of the remaining variables are held constant at their averages in the original database.

The function's gradient is computed and applied in Eq. (<ref>) and Eq. (<ref>), depending on the specific case.

Subsequently, the system is solved using a least-squares optimization method that utilizes the Trust Region Reflective [26] algorithm, which is suitable also for under-determined (non closed) systems, that may arise during the process of identifying sets of groups.

This strategy is repeated with 20 different values of $\mathbf{x}_{0}$ and iterated multiple times, each time updating the initial guess to the mean value of the exponents $(\overline{\alpha_1}, \overline{\alpha_2}, \dots, \overline{\alpha_p})$ averaged over the 20 solutions found in the previous iteration.

An overview of the methodology for identifying invariant groups/sets is depicted in Fig. <ref>.

Overview of the procedure for identifying invariant groups/sets. A regression model is trained on the physical data and used to compute the gradient of the objective function in a point $\mathbf{x}_{0}$. The matrix B is constructed according to the functional structure of the investigated group/set, and its kernel K is computed. Finally, the condition of invariance between gradient and kernel is coupled to the normalization conditions of the coefficients. If the resulting non-linear system is satisfied for the same coefficients over the $f(\mathbf{x})$ domain, the group/set is an intrinsic variable and $f(\mathbf{x})$ is invariant with respect to it.

§.§ Searching for good variables by classification models

The methodology introduced above, based on regression models, is able to identify groups and/or group sets of physical variables upon which a generic function $f(\mathbf{x})$ depends.

If successful, this approach efficiently reduces the number of such primitive variables by optimizing their combination.

However, a similar goal can be achieved by means of an entirely different approach based on classification models, and it has been initially proposed by Trezza & Chiavazzo in [23].

The latter methodology, briefly reviewed in the following, enables the optimal mixing of primitive features by performing classification and subsequent optimal class separation of the available data samples.

Let $x_{1},\dots,x_{n}$ denote $n$ features.

Let $(\Tilde{x}_{1},\dots,\Tilde{x}_{n})$ be the corresponding dimensionless quantities

\begin{equation}
\label{eq15}
\Tilde{x}_{i} = \frac{x_{i}-x_{i,\mathrm{min}}}{x_{i,\mathrm{max}}-x_{i,\mathrm{min}}}+1
\end{equation}

normalized by construction within the range $[1, 2]$ to avoid singularities in the procedure below, where $x_{i,\mathrm{min}}$ and $x_{i,\mathrm{max}}$ represent the minimum and the maximum observed values for the $i$th feature over the training set, respectively.

It is possible to synthetically create a set of $m\ll n$ mixed features $(y_{1},\dots,y_{m})$ as

\begin{equation}
\label{eq16}
y_{j}=\prod_{i=1}^{n} \Tilde{x}_{i}^{\alpha_{ij}}
\end{equation}

with $\{\alpha_{ij}\}$ being a $(n \times m)$ matrix estimated by means of a multi-objective optimization algorithm, as described below.

Finally, the new variables $y_{j}$ can be normalized within the interval $[0, 1]$ according to

\begin{equation}
\label{eq17}
\Tilde{y}_{j} = \frac{y_{j}-y_{j,\mathrm{min}}}{y_{j,\mathrm{max}}-y_{j,\mathrm{min}}}.
\end{equation}

The main idea is that the matrix $\alpha_{ij}$ shall be selected following an optimization criterion attempting the largest possible separation between two (or more) different classes partitioning the values of a physical property of interest.

Clearly, this can be accomplished by maximizing a certain distance between the classes.

However, a multi-objective optimization procedure could also be pursued.

For instance in a binary classification, the matrix $\alpha_{ij}$ in Eq. (<ref>) can be chosen to lie on the Pareto front while simultaneously pursuing: (i) maximization of a carefully selected distance between the two classes; (ii) minimization of the norm of the covariance matrix of the first class distribution; (iii) minimization of a norm of the covariance matrix of the second class distribution.

The main rationale behind the minimization of a norm of the covariance matrix for class distribution (along with a distance between classes) is the aim of possibly obtaining smooth distribution functions that can be analytically best fitted.

Following the latter idea, below we will pursue multi-objective optimization for all the examined cases and provide some optimal solutions from Pareto fronts.

In this study, we adopt genetic algorithms for optimization whereas the Bhattacharyya distance between the histograms of the two equally binned classes [27, 28] is evaluated to be maximized during the multi-objective optimization.

However, as shown in ref. [23], other options for class distance may be considered, such as the Wasserstein distance [29] or the averaged number within a fixed radius of nearest neighbors of one class to each samples of the other class.

Furthermore, herein we utilize variable power-form mixingoutlined in Eq. (<ref>). Nonetheless, alternative grouping options (e.g., linear mixing as employed in ref. [23]) are also possible with the overall methodology remaining unchanged.

Also, the procedure can be further generalized to a number of classes $>2$.

In this case, the genetic algorithm aims at simultaneously maximizing the pairwise distances between the classes [23].

For the remaining two objectives, in one-dimensional cases, a numerical estimate of the standard deviation of the binned data in the two classes is computed.

In two-dimensional (or higher) cases, the determinant of the covariance matrix can be utilized.

From the practical standpoint, considering a dataset of physical data samples characterized by features $x_1, \dots, x_n$ and their corresponding response $f(\mathbf{x})$, such dataset is partitioned into two classes based on a carefully chosen threshold for $f(\mathbf{x})$.

After the aforementioned pre-processing steps, a genetic algorithm optimization is employed to identify the Pareto front.

This optimization concurrently seeks to minimize the variances of the two classes (in one dimension) or the determinants of the covariances of the two classes (in two or more dimensions), while also maximizing the Bhattacharyya distance between the classes.

A summary of this procedure is illustrated in Fig. <ref>

Overview of the procedure to identify optimal mixed variables for class separation. Threshold values are chosen to divide the physical data in classes. A Pareto optimization is performed to construct a reduced set of synthetic features that simultaneously maximizes the Bhattacharyya distance between the classes and (i) minimizes the variances of the class distributions, in the one dimensional case, or (ii) minimizes the determinants of the covariance matrix of the class distributions, in the multi dimensional case.

§ NUMERICAL EXAMPLES AND DISCUSSION

§.§ Dittus-Boelter equation

The Dittus-Boelter correlation for a fluid undergoing heating is expressed by the following equation:

\begin{equation}
\label{eq18}
\mathrm{Nu} = 0.023\mathrm{Re}_{D}^{4/5}\mathrm{Pr}^{0.4}
\end{equation}

where $\mathrm{Re}_{D} = \frac{ud}{\nu}$ represents the Reynolds number and $\mathrm{Pr} = \frac{\nu}{\kappa}$ denotes the Prandtl number.

The equation can be rewritten by substituting the dimensionless quantities $Re$ and $Pr$, thus obtaining

\begin{equation}
\label{eq19}
\mathrm{Nu} = 0.023\left(\frac{ud}{\nu}\right)^{4/5}\left(\frac{\nu}{\kappa}\right)^{0.4} = 0.023\frac{u^{0.8}d^{0.8}}{\nu^{0.4}\kappa^{0.4}}.
\end{equation}

It is easy to verify that the objective value is invariant with respect to all the possible combinations of input features, $u, d, \nu, \kappa$, either couples or triplets.

§.§.§ Use of regression models

We attempt to recover the symmetries of the Nusselt number with respect to binary and ternary groups in the form $x_{i}^{\alpha_1}x_{j}^{\alpha_2}$ and $x_{i}^{\alpha_1}x_{j}^{\alpha_2}x_{k}^{\alpha_3}$ only relying on noised data and by adopting the methodology illustrated in subsection <ref>.

As outlined above, this requires the evaluation of the gradient of the noised Nusselt number with respect to the input features, namely $\nabla \overline{\mathrm{Nu}}(u,d,\nu,\kappa)$, where for each sample $i$, $\overline{\mathrm{Nu}}_i = \mathrm{Nu}_i(1+\eta_i)$ (see Methods for details). This can be conveniently computed by automatic differentiation over a DNN model approximating $\overline{\mathrm{Nu}}(u,d,\nu,\kappa)$.

As input features of the DNN, the four variables $u,d,\nu,\kappa$ are used.

The dataset is thus divided into three parts: (i) a training set, (ii) a validation set used to detect potential overfitting, and (iii) a testing set for the comprehensive evaluation of model performance.

Fig. <ref>a showcases the model predictions over the testing set, while in Fig. <ref>b the corresponding loss across epochs is depicted.

Notably, the model is highly predictive, with coefficient of determination $R^2=0.964$ over the testing set, with no evidence of overfitting observed.

Results of the DNN regression model for the noised Nusselt number $\overline{\mathrm{Nu}}$ in the Dittus-Boelter correlation. a Predictions over the testing set, and b corresponding loss curves for the DNN model. Model performances are shown in terms of coefficient of determination $R^{2}$, mean absolute error (MAE), and root mean squared error (RMSE).

The model is thus fed to the optimization algorithm.

For each of the ten expected invariant groups (six binary and four ternary), the algorithm estimates the normalized exponents $\alpha_1, \alpha_2, \alpha_3$ 20 times.

Table <ref> shows such true normalized exponents $\alpha_1, \alpha_2, \alpha_3$, together with the means $\mu_{\alpha_1}, \mu_{\alpha_2}, \mu_{\alpha_3}$ and the standard deviations $\sigma_{\alpha_1}, \sigma_{\alpha_2}, \sigma_{\alpha_3}$ over those 20 evaluations.

Our findings indicate that the method correctly identifies invariant groups in both couple and triplet forms.

Furthermore, it demonstrates the ability to estimate normalized exponents $\alpha_1, \alpha_2, \alpha_3$ for individual primitive variables within each group, with relative percent errors not exceeding 6.0% for couples and 16.7% for triplets.

Clearly, we notice that normalized exponents are obtained up to the sign: the exponents determined for the pair $(u, d)$ are negative instead of the expected positive values according to the Dittus-Boelter equation.

In Table <ref>, we label as found all the groups with standard deviations $\sigma_{\alpha_1}, \sigma_{\alpha_2}, \sigma_{\alpha_3} \leq 0.2$, since low values of $\sigma$ imply that the results for the 20 evaluations are consistent.

We also label as reliable all the groups where none of the evaluated exponents approach 0 or 1, since such cases would correspond to trivial solutions.

Remarkably, all six binary groups and all four ternary groups of input features are found to comply with the condition of invariance, and the results for all of them are reliable.

Normalized exponents $\alpha_1, \alpha_2, \alpha_3$ for the Dittus-Boelter correlation together with their average estimates over 20 evaluations $\mu_{\alpha_1}, \mu_{\alpha_2}, \mu_{\alpha_3}$ and the corresponding standard deviations $\sigma_{\alpha_1}, \sigma_{\alpha_2}, \sigma_{\alpha_3}$. Found groups refer to low standard deviation, reliable groups refer to average far from 1 and 0.

Group     $\alpha_1$     $\alpha_2$     $\alpha_3$     ${\mu_{\alpha_1}}$     ${\sigma_{\alpha_1}}$      ${\mu_{\alpha_2}}$     ${\sigma_{\alpha_2}}$       ${\mu_{\alpha_3}}$     ${\sigma_{\alpha_3}}$     Found     Reliable

$(u,d)$     0.71     0.71     -     -0.67     0.12     -0.73     0.10     -     -     yes     yes

$(u,\nu)$     0.89     -0.45     -     0.90     0.03     -0.44     0.07     -     -      yes     yes

$(u,\kappa)$     0.89     -0.45     -     0.90     0.02     -0.42     0.04     -     -     yes     yes

$(d,\nu)$     0.89     -0.45     -     0.88     0.06     -0.46     0.06     -     -     yes     yes

$(d,\kappa)$     0.89     -0.45     -     0.88     0.04     -0.47     0.07     -     -     yes     yes

$(\nu,\kappa)$     -0.71     -0.71     -     -0.70     0.09     -0.71     0.09     -     -     yes     yes

$(u,d,\nu)$     0.67     0.67     -0.33     0.67     0.04     0.65     0.05     -0.34     0.09     yes     yes

$(u,d,\kappa)$     0.67     0.67     -0.33     0.66     0.03     0.66     0.05     -0.34     0.07     yes     yes

$(u,\nu,\kappa)$     0.82     -0.41     -0.41     0.88     0.05     -0.34     0.09     -0.35     0.07     yes     yes

$(d,\nu,\kappa)$     0.82     -0.41     -0.41     0.86     0.03     -0.37     0.10      -0.34     0.07     yes     yes

§.§.§ Use of classification models

In a second attempt of reducing the number of input variables, we investigate the possible existence of symmetries for classification, aiming at the construction of $m$ mixed features in the form $y_{j}=\prod_{i=1}^{n} \Tilde{x}_{i}^{\alpha_{ij}}$, where $(\Tilde{x}_{1},\dots,\Tilde{x}_{n})$ are the primitive variables properly normalized within the interval [1, 2] and $\{\alpha_{ij}\}$ denotes an $n \times m$ matrix optimally estimated, as detailed in subsection <ref>.

Such features are finally properly normalized as $\Tilde{y}_j$ in the interval [0, 1].

In the case of the Dittus-Bolter, we set class 1 for samples with $\overline{\mathrm{Nu}}<395$ and class 2 for samples with $\overline{\mathrm{Nu}}\geq395$.

We thus create a single mixed feature ($m=1$) according to the above procedure.

Fig. <ref>a reports the PDF binning of the training set data over the two classes, against the normalized flow velocity $u_{\textrm{norm}}$, while Fig. <ref>b shows the same PDFs against the normalized mixed feature $\Tilde{y}_1$.

It is noteworthy to observe that when represented against the primitive feature, the two classes exhibit considerable overlap, whereas the same two classes appear well separated when plotted against the mixed feature.

As a result, it is particularly convenient to attempt an analytical best-fitting of the two distributions depicted in Fig. <ref>b approximated by a Generalized Extreme Value (GEV) distribution, with density

\begin{equation}
\label{eq21}
g(\Tilde{y}_{1}) = \frac{1}{\sigma} \left( 1+\omega \frac{\Tilde{y}_{1}-\mu}{\sigma} \right)^{-\frac{\omega+1}{\omega}}  \exp{\left( -\left( 1+\omega \frac{\Tilde{y}_{1}-\mu}{\sigma} \right)^{-1/\omega} \right)}.
\end{equation}

The fitting is performed by means of the SciPy Python package [30].

The specific computed GEV distribution for samples with $\overline{\mathrm{Nu}}<395$ has factors $\mu = 0.338$, $\sigma = 0.137$, and $\omega = 0.342$, whereas the GEV distribution for samples with $\overline{\mathrm{Nu}} \geq 395$ has factors $\mu = 0.675$, $\sigma = 0.077$, and $\omega = 0.074$.

Fig. <ref>c shows the PDFs over the binned data of the testing set reported against the same mixed feature $\Tilde{y}_{1}$, along with the GEV fittings computed on the training set.

Notably, the classes are still well separated, with a good agreement between the GEV distributions and the testing set densities.

One dimensional example for classification on the Dittus-Boelter correlation. a PDFs over binned data of the training set for the two classes ($\overline{\mathrm{Nu}}<395$ and $\overline{\mathrm{Nu}} \geq 395$) reported against the normalized flow velocity. b PDFs over binned data of the training set for the two classes reported against the mixed feature $y_{1}$, constructed according to Eq. (<ref>) and choosing the point of the Pareto front with the least overlapping of the two classes according to the Bhattacharyya distance, along with a GEV analytical fitting of the two binnings. c PDFs over binned data of the testing set for the two classes reported against the same mixed feature $y_{1}$ together with the same GEV fittings of the b subfigure. The mixed variable $\Tilde{y}_1$ shown here is referred exclusively to this Dittus-Boelter one dimensional optimization.

Furthermore, we make an attempt at the creation of two mixed normalized features $\Tilde{y}_1, \Tilde{y}_2$ ($m=2$) with the same classes.

Specifically, Fig. <ref>a shows the PDF two dimensional binning of the training set data over the two classes, against the normalized flow velocity $u_\mathrm{norm}$ and  the hydraulic diameter $d_\mathrm{norm}$.

Fig. <ref>b shows the same PDFs against the normalized mixed features $\Tilde{y}_{1}, \Tilde{y}_{2}$ constructed according to Eq. (<ref>) by power combination of the four primitive variables and choosing the point of the Pareto front with the least overlapping of the two distributions.

Similarly to the one dimensional case, the classes appear well separated when plotted against the mixed features, whereas there is a wide overlap when plotted against the primitive variables.

Fig. <ref>c shows the PDFs over the binned data of the testing set reported against the same mixed features.

The two classes appear still well separated.

Two dimensional example for the classification on the Dittus-Boelter correlation. a PDFs over binned data of the training set for the two classes ($\overline{\mathrm{Nu}}<395$ and $\overline{\mathrm{Nu}} \geq 395$) reported against the normalized flow velocity $u_{\mathrm{norm}}$ and the normalized hydraulic diameter $d_{\mathrm{norm}}$. b PDFs over binned data of the training set for the two classes reported against the mixed features $\Tilde{y}_{1}, \Tilde{y}_{2}$, constructed according to Eq. (<ref>) and choosing the point of the Pareto front with the least overlapping of the two classes according to the Bhattacharyya distance. c PDFs over binned data of the testing set for the two classes reported against the same mixed features $\Tilde{y}_{1}, \Tilde{y}_{2}$. The mixed variables $\Tilde{y}_1, \Tilde{y}_2$ shown here are referred exclusively to this Dittus-Boelter two dimensional optimization.

Finally, we aim at the construction of $m=1$ mixed feature for separating the data samples into three classes: class 0 for $\overline{\mathrm{Nu}}<197.5$, class 1 for $197.5\leq\overline{\mathrm{Nu}}<395$, class 2 for $\overline{\mathrm{Nu}}\geq395$. The multi-objective optimization is performed aiming at concurrently maximizing the pairwise distances between the classes.

Fig. <ref> shows the separation in classes reported against the normalized flow velocity $u_\mathrm{norm}$ and the new mixed feature $\Tilde{y}_{1}$: in line with previous cases, the substantial overlap observed when representing classes against the primitive feature is significantly reduced when using the optimized mixed feature.

The bet-fitting GEV distributions for the three classes are computed according to Eq. (<ref>) over the training set samples.

Specifically the GEV associated with Class 0 has factors $\mu = 0.284$, $\sigma = 0.103$, $\omega = 0.494$, the GEV for Class 1 has factors $\mu = 0.475$, $\sigma = 0.067$, $\omega = 0.193$, and the GEV for Class 2 has factors $\mu = 0.684$, $\sigma = 0.080$, $\omega = 0.103$.

When plotting the PDFs over the binned data of the testing set against the mixed feature, the classes appear well separated, and there is good agreement between the GEV distributions and the testing set densities.

One dimensional example for the ternary classification on the Dittus-Boelter. a PDFs over binned data of the training set for the three classes ($\overline{\mathrm{Nu}}<197.5$, $197.5 \leq \overline{\mathrm{Nu}} < 395$, and  $\overline{\mathrm{Nu}}\geq 395$) reported against the normalized flow velocity $u_{\mathrm{norm}}$. b PDFs over binned data of the training set for the three classes reported against the mixed feature $\Tilde{y}_{1}$, constructed according to Eq. (<ref>) and choosing the point of the Pareto front with the least overlapping of the three classes according to the Bhattacharyya distance, along with a GEV analytical fitting of the three binnings. c PDFs over binned data of the testing set for the three classes reported against the same mixed feature $\Tilde{y}_{1}$ together with the same GEV fittings of the b subfigure. The mixed variable $\Tilde{y}_1$ shown here is referred exclusively to this Dittus-Boelter one dimensional optimization.

§.§ Gnielinski correlation

The Gnielinski correlation for turbulent flow in tubes is expressed by the equation

\begin{equation}
\label{eq22}
\mathrm{Nu} = \frac{(f/8)(\mathrm{Re}_{D}-1000)\mathrm{Pr}}{1+12.7(f/8)^{1/2}(\mathrm{Pr}^{2/3}-1)}
\end{equation}

where $\mathrm{Re}_{D} = \frac{ud}{\nu}$ represents the Reynolds number, $\mathrm{Pr} = \frac{\nu}{\kappa}$ is the Prandtl number and $f$ denotes the friction factor.

The equation can be reformulated by substituting the dimensionless quantities $Re$ and $Pr$, resulting in:

\begin{equation}
\label{eq23}
\mathrm{Nu} = \frac{(f/8) \left(\frac{ud}{\nu}-1000 \right)\left(\frac{\nu}{\kappa}\right)}{1+12.7(f/8)^{1/2}\left(\left(\frac{\nu}{\kappa}\right)^{2/3}-1\right)}.
\end{equation}

The Nusselt number is invariant only with respect to the combination of flow velocity and hydraulic diameter $(u,d)$ with real exponents $(1,1)$, normalized exponents $(0.707,0.707)$.

§.§.§ Use of regression models

Here, we first attempt the detection of possible symmetries of the noised Nusselt number $\overline{\mathrm{Nu}}$ adopting the methodology proposed in subsection <ref>.

We get access to the gradient $\nabla \overline{\mathrm{Nu}}(u,d,\nu,\kappa,f)$ by means of automatic differentiation over a DNN approximating $\overline{\mathrm{Nu}}(u,d,\nu,\kappa,f)$.

Fig. <ref>a shows the model predictions over the testing set, while in Fig. <ref>b the corresponding loss across epochs is depicted.

The model is highly predictive, with coefficient of determination $R^2=0.976$ over the testing set, and no evidence of overfitting is observed.

Results of the DNN regression model for the noised Nusselt number $\overline{\mathrm{Nu}}$ in the Gnielinski correlation. a Predictions over the testing set, and b corresponding loss curves for the DNN model. Model performances are shown in terms of coefficient of determination $R^{2}$, mean absolute error (MAE), and root mean squared error (RMSE).

The trained model is thus fed to the optimization algorithm, looking for binary groups in the form $x_i^{\alpha_1} x_j^{\alpha_2}$.

The average $\mu_{\alpha_1}, \mu_{\alpha_2}$ estimates of the exponents $\alpha_1, \alpha_2$, together with their standard deviations  $\sigma_{\alpha_1}, \sigma_{\alpha_2}$ are reported in Table <ref>.

Specifically, the procedure correctly identifies the expected group $(u,d)$, whereas for most other groups the results correspond to either the trivial solution (i.e., the absolute value of one of the exponents is close to 1 and the other is close to 0), or at least one of the evaluated values presents too high variance (i.e., $\sigma_{\alpha_i} > 0.2$), as for $(\kappa,f)$.

Remarkably, for the remaining cases (namely $(u,\nu)$, $(u,f)$, $(d,\nu)$, $(d,f)$, $(\nu,\kappa)$, $(\nu,f)$) the procedure identifies groups complying with the invariance condition, albeit not explicitly expected in the Gnielinski correlation.

We thus concluded that the suggested procedure has detected a local invariance and this can be numerically verified as follows (here we limit to groups of two variables only).

First, the features not appearing in the group are kept constant.

Second, a random sample $i$ in the dataset is considered for the evaluation of the quantity $\Tilde{c}=x_{1,i}^{\alpha_1} x_{2,i}^{\alpha_2}$.

Third, a vector $\overline{\mathbf{x}_{1}}$ of evenly spaced points in the domain of the first feature is created.

As such, the array corresponding to the second feature is obtained as $\overline{\mathbf{x}_2} = (\Tilde{c}\overline{\mathbf{x}_1}^{-\alpha_1})^{1/\alpha_2}$.

A new dataset is thus constructed, with the variables out of the group being constant and with the variables in the group being replaced by $\overline{\mathbf{x}_1}$, $\overline{\mathbf{x}_2}$.

For all those samples, the response value $f(\mathbf{x})$ is computed.

The local invariance is demonstrated when $f(\mathbf{x})$ is approximately constant over all the new constructed dataset.

Further details can be found in Supplementary Note 3.

89Normalized exponents $\alpha_1, \alpha_2, \alpha_3, \beta_1, \beta_2$ for the Gnielinski correlation, together with their average estimates over 20 evaluations $\mu_{\alpha_1}$, $\mu_{\alpha_2}$, $\mu_{\alpha_3}$, $\mu_{\beta_1}$, $\mu_{\beta_2}$ and the corresponding standard deviations $\sigma_{\alpha_1}$, $\sigma_{\alpha_2}$, $\sigma_{\alpha_3}$, $\sigma_{\beta_1}$, $\sigma_{\beta_2}$. Found groups/sets refer to low standard deviation, reliable groups refer to average far from 1 and 0.

Group/Set     $\alpha_1$     $\alpha_2$     $\alpha_3$     $\beta_1$     $\beta_2$     ${\mu_{\alpha_1}}$     ${\sigma_{\alpha_1}}$     ${\mu_{\alpha_2}}$     ${\sigma_{\alpha_2}}$     ${\mu_{\alpha_3}}$     ${\sigma_{\alpha_3}}$     ${\mu_{\beta_1}}$     ${\sigma_{\beta_1}}$     ${\mu_{\beta_2}}$     ${\sigma_{\beta_2}}$     Found     Reliable

$(u,d)$     0.707     0.707     -     -     -     -0.727     0.050     -0.682     0.055     -     -     -     -     -     -     yes     yes

$(u,\nu)$     -     -     -     -     -     0.830     0.041     -0.553     0.062     -     -     -     -     -     -     yes     yes

$(u,\kappa)$     -     -     -     -     -     0.944     0.025     -0.322     0.072     -     -     -     -     -     -     yes     no

$(u,f)$     -     -     -     -     -     -0.829     0.082     -0.521     0.132     -     -     -     -     -     -     yes     yes

$(d,\nu)$     -     -     -     -     -     0.838     0.045     -0.538     0.073     -     -     -     -     -     -     yes     yes

$(d,\kappa)$     -     -     -     -     -     -0.928     0.053     0.350     0.117     -     -     -     -     -     -     yes     no

$(d,f)$     -     -     -     -     -     0.862     0.080     0.457     0.104     -     -     -     -     -     -     yes     yes

$(\nu,\kappa)$     -     -     -     -     -     -0.896     0.065     -0.405     0.170     -     -     -     -     -     -     yes     yes

$(\nu,f)$     -     -     -     -     -     0.720     0.070     -0.687     0.075     -     -     -     -     -     -     yes     yes

$(\kappa,f)$     -     -     -     -     -     0.518     0.142     -0.839     0.285     -     -     -     -     -     -     no     -

$[(u,\nu),(\nu,\kappa)]$     0.707     -0.707     -     0.707     -0.707     0.709     0.066     -0.698     0.067     -     -     -0.692     0.017     0.722     0.016     yes     yes

$[(d,\nu),(\nu,\kappa)]$     0.707     -0.707     -     0.707     -0.707     0.702     0.073     -0.704     0.076     -     -     -0.700     0.019     0.714     0.019     yes     yes

$[(u,d),(u,\nu)]$     -     -     -     -     -     0.260     0.042     -0.965     0.010     -     -     -0.903     0.055     0.314     0.287     yes     no

$[(u,d),(d,f)]$     -     -     -     -     -     -0.541     0.059     -0.838     0.042     -     -     -0.569     0.020     0.822     0.014     yes     yes

$[(u,\nu),(\nu,f)]$     -     -     -     -     -     0.858     0.024     -0.512     0.044     -     -     -0.183     0.047     0.982     0.006     yes     no

$[(u,f),(f,d)]$     -     -     -     -     -     -0.632     0.080     -0.768     0.069     -     -     -0.495     0.058     0.867     0.035     yes     yes

$[(d,f),(f,\nu)]$     -     -     -     -     -     0.849     0.018     -0.527     0.028     -     -     -0.883     0.038     0.463     0.064     yes     yes

$[(u,d,\nu),(\nu,\kappa)]$     0.577     0.577     -0.577     0.707     -0.707     0.587     0.045     0.562     0.073     -0.572     0.078     0.700     0.023     -0.713     0.022     yes     yes

$[(u,d,f),(f,\nu)]$     -     -     -     -     -     0.573     0.005     0.544     0.031     -0.612     0.032     0.833     0.023     -0.551     0.033     yes     yes

$[(u,\kappa,f),(\kappa, \nu)]$     -     -     -     -     -     0.422     0.055     0.721     0.034     -0.540     0.077     0.606     0.029     -0.794     0.021     yes     yes

Furthermore, armed with the same DNN regression model, we try to detect invariance of the Nusselt number with respect to sets groups in the form $x_{i}^{\alpha_1}x_{j}^{\alpha_2}x_{k}^{\alpha_3}$ and $x_{k}^{\beta_1}x_{l}^{\beta_2}$ employing the procedure presented in subsection <ref>.

<Ref> reports also the results for selected sets of two feature couples, including sets appearing explicitly in the Gnielinski correlation – i.e., $[(u,\nu),(\nu,\kappa)]$ and $[(d,\nu),(\nu,\kappa)]$ – and for selected sets of features comprising a triplet and a couple – e.g., $[(u,d,\nu),(\nu,\kappa)]$, which is extremely relevant as it comprises the Reynolds and the Prandtl numbers.

For all the mentioned sets the procedure correctly identifies the exponents that make the groups compliant with the condition of invariance.

Indeed, the normalized exponents obtained with the procedure for the sets $[(u,\nu),(\nu,\kappa)]$ and $[(d,\nu),(\nu,\kappa)]$ are $[(0.709, -0.698), (-0.692, 0.722)]$ and $[(0.702, -0.704), (-0.700, 0.714)]$ respectively, whereas the analytical solution is $[(0.707,-0.707),(0.707,-0.707)]$ for both the cases.

The normalized exponents obtained for the set $[(u,d,\nu),(\nu,\kappa)]$ are $[(0.587, 0.562, -0.572), (0.700, -0.713)]$, whereas the analytical solution is $[(0.577,0.577,-0.577),(0.707,-0.707)]$.

<Ref> also reports some identified sets that locally comply with the condition of invariance (see Supplementary Note 3 for more details), together with some randomly selected sets for which the procedure does succeed at finding the invariance property, but for which the results are close to the trivial solution, meaning that at least one of the evaluated exponents approaches either 0 or $\pm1$, thus being not reliable.

Following the same approach above, Table <ref> labels as found groups with low standard deviations and as reliable all the groups in which none of the evaluated exponents approaches 0 or $\pm1$.

Remarkably, the analytically evident couple $(u, d)$, sets of two couples $[(u,\nu),(\nu,\kappa)]$ and $[(d,\nu),(\nu,\kappa)]$, and set comprising a triplet and a couple $[(u,d,\nu),(\nu,\kappa)]$, are identified.

Furthermore, eight additional couples are found that locally comply with the condition of invariance, six of them being reliable, and only one couple is not found.

In the case of the selected sets of two feature couples, five more sets of two couples locally comply with the condition of invariance, only two of them not being reliable.

Finally, the two additional sets comprising a triplet and a couple are reliably found to represent a local invariance.

§.§.§ Use of classification models

In a second attempt at reducing the number of input variables, we aim at the construction of mixed optimized features as power combinations of the normalized primitive variables for classification, according to Eq. (<ref>).

In the case of Gnielinski correlation, we set class 1 for samples with $\overline{\mathrm{Nu}}<500$ and class 2 for samples with $\overline{\mathrm{Nu}}\geq500$.

Once the optimization routine finds the Pareto front, the mixed features are created using the point of the Pareto front with the least overlapping of the two classes according to the Battacharyya distance.

Fig. <ref>a reports the PDF binning of the training set data over the two classes, against the normalized flow velocity $u_{\textrm{norm}}$, while Fig. <ref>b shows the same PDFs against the normalized mixed feature $\Tilde{y}_1$; still, when represented against the primitive variable, the two classes exhibit significant overlap, whereas they appear distinctly separated when plotted against the mixed feature.

Two bet-fitting GEV distributions are computed for the two classes according to Eq. (<ref>) for samples in class 1 with factors $\mu = 0.276$, $\sigma = 0.112$, $\omega = 0.223$; for samples in class 2 with factors $\mu = 0.510$, $\sigma = 0.091$, $\omega = 0.000$. Fig. <ref>c shows the PDFs over the binned data of the testing set reported against the same mixed feature $\Tilde{y}_{1}$, along with the GEV fittings computed on the training set. Notably, the classes are still well separated, with a good agreement between the GEV distributions and the testing set densities.

One dimensional example: a PDFs over binned data of the training set for the two classes ($\overline{\mathrm{Nu}}<500$ and $\overline{\mathrm{Nu}}\geq 500$) reported against the normalized flow velocity. b PDFs over binned data of the training set for the two classes reported against the mixed feature $y_{1}$, constructed according to Eq. (<ref>) and choosing the point of the Pareto front with the least overlapping of the two classes according to the Bhattacharyya distance, along with a GEV analytical fitting of the two binnings. c PDFs over binned data of the testing set for the two classes reported against the same mixed feature $y_{1}$ together with the same GEV fittings of the b subfigure. The mixed variable $\Tilde{y}_1$ shown here is referred exclusively to this Gnielinski one dimensional optimization.

Furthermore, we attempt to create two mixed normalized features $\Tilde{y}_1, \Tilde{y}_2$ with the same classes. Fig. <ref>a shows the PDF two dimensional binning of the training set data over the two classes, against the normalized flow velocity $u_\mathrm{norm}$ and friction factor $f_\mathrm{norm}$. Fig. <ref>b shows the same PDFs against the mixed features $\Tilde{y}_{1}, \Tilde{y}_{2}$ constructed according to Eq. (<ref>) by power combination of the five relevant features and choosing the point of the Pareto front with the least overlapping of the classes. Similarly to the one dimensional case, the classes appear again well separated when plotted against the mixed features, whereas there is a wide overlap when plotted against the primitive variables. <Ref>c shows the PDFs over the binned data of the testing set, reported against the same mixed features $\Tilde{y}_{1}, \Tilde{y}_{2}$; the two classes still appear well separated.

Two dimensional example for the classification on the Gnielinski correlation: a PDFs over binned data of the training set for the two classes ($\overline{\mathrm{Nu}}<500$ and $\overline{\mathrm{Nu}}\geq 500$) reported against the normalized flow velocity $u_{\mathrm{norm}}$ and friction factor $f$. b PDFs over binned data of the training set for the two classes reported against the mixed features $\Tilde{y}_{1}, \Tilde{y}_{2}$, constructed according to Eq. (<ref>) and choosing the point of the Pareto front with the least overlapping of the two classes according to the Bhattacharyya distance. c PDFs over binned data of the testing set for the two classes reported against the same mixed features $\Tilde{y}_{1}, \Tilde{y}_{2}$. The mixed variables $\Tilde{y}_1, \Tilde{y}_2$ shown here are referred exclusively to this Gnielinski two dimensional optimization.

Finally, we aim at the construction of $m=1$ mixed feature for separating the data samples into three classes: class 0 for $\overline{\mathrm{Nu}}<400$, class 1 for $400\leq\overline{\mathrm{Nu}}<900$, class 2 for $\overline{\mathrm{Nu}}\geq 900$. The multi-objective optimization is performed aiming at concurrently maximizing the pairwise distances between the classes.

Fig. <ref> shows the separation in classes reported against the normalized flow velocity $u_\mathrm{norm}$ and the new mixed feature $\Tilde{y}_{1}$: in line with previous cases, the considerable overlap observed when representing classes against the primitive feature is significantly reduced when using the optimized mixed feature.

The bet-fitting GEV distributions for the three classes are computed according to Eq. (<ref>) over the training set samples; specifically the GEV associated with Class 0 has factors $\mu = 0.245$, $\sigma = 0.083$, $\omega = 0.346$, the GEV for Class 1 has factors $\mu = 0.430$, $\sigma = 0.053$, $\omega = 0.143$, and the GEV for Class 2 has factors $\mu = 0.599$, $\sigma = 0.074$, $\omega = 0.017$. When plotting the PDFs over the binned data of the testing set against the mixed feature, the classes appear well separated, and there is good agreement between the GEV distributions and the testing set densities.

One dimensional example for the ternary classification on the Gnielinski correlation. a PDFs over binned data of the training set for the three classes ($\overline{\mathrm{Nu}}<400$, $400 \leq \overline{\mathrm{Nu}} < 900$, and  $\overline{\mathrm{Nu}}\geq 900$) reported against the normalized flow velocity $u_{\mathrm{norm}}$. b PDFs over binned data of the training set for the three classes reported against the mixed feature $\Tilde{y}_{1}$, constructed according to Eq. (<ref>) and choosing the point of the Pareto front with the least overlapping of the three classes according to the Bhattacharyya distance, along with a GEV analytical fitting of the three binnings. c PDFs over binned data of the testing set for the three classes reported against the same mixed feature $\Tilde{y}_{1}$ together with the same GEV fittings of the b subfigure. The mixed variable $\Tilde{y}_1$ shown here is referred exclusively to this Gnielinski one dimensional optimization.

§.§ Newton's law of universal gravitation

The module of the interacting force $F_\mathrm{g}$ for two objects of masses $m_1, m_2$ is expressed via Newton's law of universal gravitation

\begin{equation}
\label{eq}
F_{\mathrm{g}} = \frac{Gm_{1}m_{2}}{(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}+(z_{1}-z_{2})^{2}},
\end{equation}

where $G$ is the gravitational constant, $(x_1,y_1,z_1)$ and $(x_2,y_2,z_2)$ are the coordinates of the centers of their masses.

In this case, our aim is to identify the invariance of the noised $\overline{F_{\mathrm{g}}}$ with respect to groups $(x_1 - x_2)$, $(y_1 - y_2)$, $(z_1 - z_2)$. To this end, we employ the procedure presented in subsection <ref> for identification of general group form, focusing specifically on groups with functional dependence $\alpha_1 x_i+\alpha_2 x_j$.

As usual, we get access to the gradient $\nabla \overline{F_\mathrm{g}}(x_1, x_2, y_1, y_2, z_1, z_2, m_1, m_2)$ by means of automatic differentiation over a DNN approximating $\overline{F_\mathrm{g}}(x_1, x_2, y_1, y_2, z_1, z_2, m_1, m_2)$.

Fig. <ref>a shows the model predictions over the testing set, while in Fig. <ref>b the corresponding loss across epochs is depicted; specifically, the model is highly predictive, with coefficient of determination $R^2=0.959$ and no evidence of overfitting is observed.

Results of the DNN regression model for the noised gravitational force $F_\mathrm{g}$. a Predictions over the testing set, and b corresponding loss curves for the DNN model. Model performances are shown in terms of coefficient of determination $R^{2}$, mean absolute error (MAE), and root mean squared error (RMSE).

The trained model is thus fed to the optimization algorithm. The average $\mu_{\alpha_1}, \mu_{\alpha_2}$ estimates of the coefficients $\alpha_1, \alpha_2$, together with their standard deviations $\sigma_{\alpha_1}, \sigma_{\alpha_2}$ are reported in <Ref>. Specifically, the procedure correctly identifies the groups $(x_1, x_2)$, $(y_1, y_2)$, $(z_1, z_2)$ with estimates $\mu_{\alpha_1}, \mu_{\alpha_2}$ of $(-0.729, 0.680)$, $(-0.715, 0.681)$, $(-0.702, 0.697)$ respectively, compared to the true normalized coefficients $(-0.707, 0.707)$ for all cases.

Table <ref> flags as found groups with low standard deviations and as reliable all the groups in which none of the evaluated exponents approaches 0 or 1. Specifically, other variables couples show high variance over the corresponding average estimates, meaning that no further invariant group is identified.

Normalized exponents $\alpha_1, \alpha_2$ for the Newton's law of universal gravitation, together with their average estimates over 20 evaluations $\mu_{\alpha_1}$, $\mu_{\alpha_2}$, and the corresponding standard deviations $\sigma_{\alpha_1}$, $\sigma_{\alpha_2}$. Found groups/sets refer to low standard deviation, reliable groups refer to average far from 1 and 0.

Group     $\alpha_1$     $\alpha_2$     ${\mu_{\alpha_1}}$     ${\sigma_{\alpha_1}}$     ${\mu_{\alpha_2}}$     ${\sigma_{\alpha_2}}$     Found     Reliable

$(x_1, x_2)$     0.707     -0.707     -0.729     0.047     0.680     0.057     yes     yes

$(y_1, y_2)$     0.707     -0.707     -0.715     0.09     0.681     0.128     yes     yes

$(z_1, z_2)$     0.707     -0.707     -0.702     0.103     0.697     0.104     yes     yes

$(m_1, m_2)$     -     -     -0.455     0.601     0.294     0.587     no     -

$(m_1, x_1)$     -     -     0.000     0.000     1.000     0.000     no     -

$(x_1, y_1)$     -     -     -0.628     0.271     -0.032     0.714     no     -

§ CONCLUSIONS AND FINAL REMARKS

In this study, we have implemented two innovative methodologies for searching optimal variables to describe physical data, making use of both regression and classification models applied to the data.

Specifically, leveraged on well-suited datasets for machine learning, with the goal of predicting a property of interest.

More precisely, the procedure based on the regression model introduced here enables the identification of invariant groups and/or sets of variable groups with respect to which the property of interest is invariant.

We have demonstrated its effectiveness on noised data generated by three well known functional relationships: the Dittus-Boelter correlation, the Gnielinski correlation, and Newton's law of universal gravitation.

For the former two cases, the procedure has accurately detected groups/sets in power form, while for the latter case, a generalized algorithm successfully identified groups with a linear form.

Interestingly, the methodology is potentially applicable to any functional form, as illustrated in subsection <ref>.

The procedure based on classification of the physical property values, allows the determination of an appropriate set of exponents to combine all the primitive variables in power form, thus constructing mixed features optimized for the classification task.

Specifically, we have shown its effectiveness on the Dittus-Boelter and on the Gnielinski correlations.

Indeed, the methodology effectively enables the separation of classes even with just one mixed feature, whereas a single primitive variable fails to achieve class separation.

Furthermore, we also provide examples with one and two mixed variables, together with separation in two or three classes.

It is worth stressing that the methodologies presented in this study are blessed by generality and, as such, are not restricted to the selected case studies.

Therefore, potential applications can be envisioned in various other fields in the future.

In particular, we notice that the identification of effective good variables is not only interesting per se to possibly gain a deeper insight on a physical system, but can also be practically advantageous when designing experiments.

Indeed, the methodology to detect groups/sets in regression facilitates efficient group/set-level adjustments rather than individually fine-tuning variables within the groups/sets themselves.

Moreover, the classification procedure can enable the reduction of the numerous original primitive variables to a minimal set of optimized variables concerning a specific physical property of interest.

Notably, combinations of variables yielding identical mixed feature values exhibit similar performance in terms of the property being classified.

As a result, it is possible to find alternative combinations of primitive variables without compromising the overall performance of a given system under study.

From a more practical standpoint, the methodologies described here can help generalizing optimal system conditions, thus helping decreasing the most resource-intensive components while properly re-balancing the others.

As an example, these general methodologies may thus hold the potential to save resources, such as costly reagents (e.g., Bonke et al. [24] for solar fuel production) or expensive materials as in perovskite solar cells optimization [31].

An additional advantage associated to the correct identification of a reduced set of ruling variables is their possible use for driving sequential learning or Bayesian optimization processes [32, 33].

Finally, we believe that an interesting development of the presented methodologies to be pursued in the future, shall be in the direction of a possible handling of systems ruled not only by numerical parameters but also categorical ones.

§ DATA AVAILABILITY

Processed datasets and trained models of this study are publicly available in Zenodo at DOI: 10.5281/zenodo.10406490

§ CODE AVAILABILITY

The codes used to obtain the results of this study are publicly available in github at .

[1]
authorRappaz, M., authorBellet, M., authorDeville, M. O. & authorSnyder, R.
titleNumerical modeling in materials science and engineering (publisherSpringer, year2003).

[2]
authorChen, B. et al.
titleAutomated discovery of fundamental variables hidden in experimental data.
journalNature Computational Science volume2, pages433–442 (year2022).

[3]
authorFloryan, D. & authorGraham, M. D.
titleData-driven discovery of intrinsic dynamics.
journalNature Machine Intelligence volume4, pages1113–1120 (year2022).

[4]
authorEva, B., authorRied, K., authorMüller, T. & authorBriegel, H. J.
titleHow a minimal learning agent can infer the existence of unobserved variables in a complex environment.
journalMinds and Machines volume33, pages185–219 (year2023).

[5]
authorChiavazzo, E.
titleApproximation of slow and fast dynamics in multiscale dynamical systems by the linearized relaxation redistribution method.
journalJournal of Computational Physics volume231, pages1751–1765 (year2012).

[6]
authorChiavazzo, E. & authorKarlin, I. V.
titleQuasi-equilibrium grid algorithm: Geometric construction for model reduction.
journalJournal of Computational Physics volume227, pages5535–5560 (year2008).

[7]
authorRayleigh, L.
titleViii. on the question of the stability of the flow of fluids.
journalThe London, Edinburgh, and Dublin philosophical magazine and journal of science volume34, pages59–70 (year1892).

[8]
authorBuckingham, E.
titleOn physically similar systems; illustrations of the use of dimensional equations.
journalPhysical review volume4, pages345 (year1914).

[9]
authorCurtis, W., authorLogan, J. D. & authorParker, W.
titleDimensional analysis and the pi theorem.
journalLinear Algebra and its Applications volume47, pages117–126 (year1982).

[10]
authorChiavazzo, E. et al.
titleIntrinsic map dynamics exploration for uncharted effective free-energy landscapes.
journalProceedings of the National Academy of Sciences volume114, pagesE5494–E5503 (year2017).

[11]
authorChiavazzo, E., authorGear, C. W., authorDsilva, C. J., authorRabin, N. & authorKevrekidis, I. G.
titleReduced models in chemical kinetics via nonlinear data-mining.
journalProcesses volume2, pages112–140 (year2014).

[12]
authorLin, K. K. & authorLu, F.
titleData-driven model reduction, wiener projections, and the koopman-mori-zwanzig formalism.
journalJournal of Computational Physics volume424, pages109864 (year2021).

[13]
authorMcRee, R. K.
titleSymbolic regression using nearest neighbor indexing.
In booktitleProceedings of the 12th annual conference companion on Genetic and evolutionary computation, pages1983–1990 (year2010).

[14]
authorStijven, S., authorMinnebo, W. & authorVladislavleva, K.
titleSeparating the wheat from the chaff: on feature selection and feature importance in regression random forests and symbolic regression.
In booktitleProceedings of the 13th annual conference companion on Genetic and evolutionary computation, pages623–630 (year2011).

[15]
authorMcConaghy, T.
titleFfx: Fast, scalable, deterministic symbolic regression technology.
journalGenetic Programming Theory and Practice IX pages235–260 (year2011).

[16]
authorArnaldo, I., authorO'Reilly, U.-M. & authorVeeramachaneni, K.
titleBuilding predictive models via feature synthesis.
In booktitleProceedings of the 2015 annual conference on genetic and evolutionary computation, pages983–990 (year2015).

[17]
authorBrunton, S. L., authorProctor, J. L. & authorKutz, J. N.
titleDiscovering governing equations from data by sparse identification of nonlinear dynamical systems.
journalProceedings of the national academy of sciences volume113, pages3932–3937 (year2016).

[18]
authorQuade, M., authorAbel, M., authorNathan Kutz, J. & authorBrunton, S. L.
titleSparse identification of nonlinear dynamics for rapid model recovery.
journalChaos: An Interdisciplinary Journal of Nonlinear Science volume28 (year2018).

[19]
authorSearson, D. P., authorLeahy, D. E. & authorWillis, M. J.
titleGptips: an open source genetic programming toolbox for multigene symbolic regression.
In booktitleProceedings of the International multiconference of engineers and computer scientists, vol. volume1, pages77–80 (organizationCiteseer, year2010).

[20]
authorDubčáková, R.
titleEureqa: software review (year2011).

[21]
authorSchmidt, M. & authorLipson, H.
titleDistilling free-form natural laws from experimental data.
journalscience volume324, pages81–85 (year2009).

[22]
authorUdrescu, S.-M. & authorTegmark, M.
titleAi feynman: A physics-inspired method for symbolic regression.
journalScience Advances volume6, pageseaay2631 (year2020).

[23]
authorTrezza, G. & authorChiavazzo, E.
titleLeveraging composition-based energy material descriptors for machine learning models.
journalMaterials Today Communications volume36, pages106579 (year2023).

[24]
authorBonke, S. et al.
titleMulti-variable multi-metric optimisation of self-assembled photocatalytic co2 reduction performance using machine learning algorithms  (year2023).

[25]
authorLide, D. R. & authorKehiaian, H. V.
titleCRC handbook of thermophysical and thermochemical data (publisherCRC press, year2020).

[26]
authorBranch, M. A., authorColeman, T. F. & authorLi, Y.
titleA subspace, interior, and conjugate gradient method for large-scale bound-constrained minimization problems.
journalSIAM Journal on Scientific Computing volume21, pages1–23 (year1999).

[27]
authorBhattacharyya, A.
titleOn a measure of divergence between two statistical populations defined by their probability distribution.
journalBulletin of the Calcutta Mathematical Society volume35, pages99–110 (year1943).

[28]
authorBhattacharyya, A.
titleOn a measure of divergence between two multinomial populations.
journalSankhyā: the indian journal of statistics pages401–406 (year1946).

[29]
authorVillani, C. et al.
titleOptimal transport: old and new, vol. volume338 (publisherSpringer, year2009).

[30]
authorVirtanen, P. et al.
titleSciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.
journalNature Methods volume17, pages261–272 (year2020).

[31]
authorTailor, N. K. et al.
titleRecent progress in morphology optimization in perovskite solar cell.
journalJournal of Materials Chemistry A volume8, pages21356–21386 (year2020).

[32]
authorHuan, X. & authorMarzouk, Y. M.
titleSimulation-based optimal bayesian experimental design for nonlinear systems.
journalJournal of Computational Physics volume232, pages288–317 (year2013).

[33]
authorMotoyama, Y. et al.
titleBayesian optimization package: Physbo.
journalComputer Physics Communications volume278, pages108405 (year2022).