LCEN: A Novel Feature Selection Algorithm for Nonlinear, Interpretable Machine Learning Models

Statistical models are powerful tools to explain, predict, or describe natural phenomena [1]. They connect independent variables (also called “inputs” or “features”) to dependent variables (also called “outputs” or “labels”) to test causal hypotheses, predict novel outputs from known inputs, or summarizing the data structure [1]. Many model architectures exist, including linear, ensemble-based, and deep learning models. Complex architectures are claimed to have greater capability to model phenomena due to their lower bias, but their intricate and numerous mathematical transformations prevent humans from understanding how an output was predicted by a model, or the relative or absolute importance of the inputs. Moreover, a lack of transparency may prevent the model from being trusted in critical or sensitive applications [2].

As summarized by Ref. [3], there are two main methods to increase interpretability: the use of model-agnostic algorithms, which extract interpretable explanations a posteriori and work for any architecture, or the direct use of interpretable architectures. Interpretable architectures include “decision trees, rules, additive models, attention-based networks, and sparse linear models” [3]. It should be noted that nonlinear models may also be made sparse, and even interpretable, as described later in this section and the rest of this work. As elaborated in the review of Ref. [4], interpretable architectures can have many advantages over black-box or a posteriori explanations, including the ability to assist researchers in refining the model and data, or better highlighting scenarios in which the model fails or lacks robustness. Special attention should be given to sparse models, which identify the most important features, can make the model more robust to variations in the input data, and can significantly improve the model’s interpretability if an interpretable architecture is used [4]. At the same time, even a linear model or decision tree/rules can become unwieldy and uninterpretable if hundreds or thousands of coefficients or rules are present.

Feature selection is the process of selecting the most important features in a model to increase its robustness or interpretability. Many criteria for feature selection exist [5], including significance based on p-values (using a univariate, iterative/stepwise, or global method), using information criteria (such as the AIC [6] and BIC [7]), using penalties (such as in LASSO [8, 9] and elastic net [EN] [10]), criteria based on changes in estimates, and expert knowledge. More broadly, these methods can be classified as filter, wrapper, or integrated methods. While no method is superior for all problems, different works have evaluated and criticized these criteria. For example, stepwise regression is one of the most commonly used methods in many fields thanks to its computational simplicity and ease of understanding [11, 5, 12]. However, stepwise regression is prone to ignoring features with causal effects, including irrelevant features, generating excessively small confidence intervals, and producing incorrect/biased parameters [11, 12]. LASSO is simple and computationally cheap, and has performed well for some problems [13, 14, 15], but can overselect irrelevant variables, can select only as many features as there are samples, and does not handle multicollinear data well [5, 10].¹¹1This last point is somewhat controversial in the literature, see Refs. [13, 16], for example.

Originally, most feature selection methods applied only in linear contexts (or have been applied primarily in linear contexts). For example, the only sparse models referenced in the highly cited review by Ref. [3] are linear. The main sparse architectures (LASSO, EN, and their variants) are linear regressors.²²2Ref. [17] has claimed that, if thousands of samples are available, LASSO can consistently select features even in nonlinear settings, although the coefficients may be distorted. However, this consistency is disputed by Ref. [18]. Later works consider sparse nonlinear models, to address the fact that many natural and industrial processes are nonlinear. Refs. [19, 20, 21] defined sets of features consisting of polynomials (all works), interactions (all works), and/or non-polynomials [19, 21]. ALVEN, the model architecture from Ref. [21], uses an F-test for each feature (including the expanded set of features) to determine whether to keep a feature in the final EN model, a filter approach. However, this F-test has very poor feature selectivity, as nearly all features are selected when traditional values of $\alpha$ ($0.001\leq\alpha\leq 0.05$) are used. Furthermore, the ordering of the features with respect to their p-values does not follow their relevance, as many irrelevant features are among those with the lowest p-values, and relevant features can be among those with the highest p-values (including p $\gg 0.05$). The work of Ref. [22] uses two forms of modified LASSO algorithms for nonlinear feature selection, achieving high sparsity with the datasets tested.

More recently, the L₁ regularization has been applied to neural networks for nonlinear feature selection. In its simplest form, group LASSO is applied to zero all the outputs of some neurons (sparsifying the network and eliminating features when zeroing input-layer neurons) [18, 23, 24]. A slightly modified version of this algorithm applies the L₁ penalty only to the input layer and includes a skip-connection between that layer and the output layer [25]. More complex applications of this method include the multi-modal neural networks of Ref. [26], the concrete autoencoders of Ref. [27], and the teacher-student network of Ref. [28]. The first two are notable for being at least partially unsupervised methods, suggesting they can select the most relevant features for a given dataset no matter the task. Some works have also used approaches other than the L₁ norm, such as the L₀ norm [29]. While these models are powerful tools for feature selection, two considerable limitations are present: first, they do not provide any information on how the selected features are contributing to the final prediction, significantly limiting interpretability. A posteriori methods to extract this information have been found unreliable [4], even if useful. Second, these complex model architectures may take “shortcuts” to make apparently accurate predictions [30, 31]. However, these “shortcuts” are not really relevant to the task, preventing proper generalization and human interpretation.

To create nonlinear, interpretable machine learning models with high predictive and descriptive power, we propose the LASSO-Clip-EN (LCEN) algorithm. This algorithm generates an expanded set of nonlinear features (such as in ALVEN) and performs feature selection and model fitting. This feature set expansion, together with the Clip step, provide LCEN with the ability to do nonlinear predictions. The algorithm is tested on artificial and empirical data, successfully rediscovering physical laws from data belonging to multiple different areas of knowledge with errors < 2% on the coefficients, a value within the empirical noise of the datasets. On datasets from processes whose underlying physical laws are not yet known, LCEN attains lower RMSEs than many sparse and dense methods, leads to sparser models than any other architecture tested, and simultaneously runs faster than most alternative architectures.

This work introduces the LASSO-Clip-EN (LCEN) algorithm for the creation of nonlinear, interpretable machine learning models (Algorithm 1). LCEN is first validated using artificial data (Section 3.1), which provide an initial assessment of the algorithm’s performance under multiple, independently controllable conditions. LCEN displayed high performance for data that are noisy (“Linear, 5-variable polynomial” and “4th-degree, univariate polynomial” datasets; Table 1 and Fig. 1), multicollinear (“Multicollinear data” dataset; Fig. A1), or generated from a process with a complex nonlinear function (“Relativistic energy” and “4th-degree, univariate polynomial” datasets; Table 2 and Fig. 1), which also leads to high hyperparameter and feature variance. For the “4th-degree, univariate polynomial” dataset, LCEN automatically adjusted its predictions to compensate for the bias-variance tradeoff (Figs. 1, A2, and A3). During the experiments with the “Relativistic energy” dataset, the range of the data affected the prediction performance – data with higher ranges led to more accurate predictions even at high noises. We hypothesize that this phenomenon occurs because it is difficult to distinguish between multiple powers of a variable (or, in a more generalized manner, how Taylor approximations are successful) when the data cover low ranges and are noisy.

LCEN was then tested with data from processes with known physical laws (Section A6.1). These data have real properties (such as noise) and may have limitations common to many processes (such as low amounts of data). Nevertheless, LCEN successfully selected the correct features for all datasets used in this work, effectively rediscovering physical laws solely from data (Table 3). The “CARMENES star data” dataset contains hundreds of points described by the Stefan-Boltzmann law. This law uses a 6th-order interaction, but since this would not be known a priori, LCEN was tested using all degree hyperparameters from 1 to 10. Despite the high variance of the process and the many potential features that could be selected, LCEN selected only the correct $R^{2}T^{4}$ feature with a coefficient relative error of only 1.74%, which is well within the 2–3% data error. Next, LCEN was evaluated for the “Kepler’s 3rd Law” datasets, which contain small amounts of data (only 6 or 8 points) yet were enough for Kepler to derive the eponymous law in 1619. LCEN can automatically replicate this discovery despite the slightly inaccurate data from Kepler’s measurements, selecting only the correct $a^{3/2}$ feature. The coefficient relative errors were 0.46% on Kepler’s original data and 0.07% on modern data, extremely small values that highlight the algorithm’s accuracy.

Lastly, experiments using empirical data from processes with unknown laws further validated LCEN’s feature selection and predictive performance. LCEN was compared to multiple different model architectures, including linear models, ensemble methods, and sparse and dense MLPs. On the “Diesel Freezing Point” dataset (Table 4), LCEN had a test RMSE equal to the RMSEs for RR and EN. However, LCEN required 11.1- and 8.3-fold fewer features than RR and EN (respectively) to reach the same accuracy, and LCEN was 3.0-fold faster than EN. LCEN also generated a more accurate and sparser model than LASSO, and LCEN models with increased sparsities typically displayed only small increases in the prediction error. When compared to SVM, the architecture with the lowest test RMSE, LCEN displayed only a 10.0% higher error, despite being a completely interpretable architecture and using only 36/401 features. LCEN provided a combination of accuracy, sparsity, interpretability, and speed that was unmatched by the other machine-learning methods. On the “Abalone” dataset (Table A7), LCEN automatically recognized the importance of 2nd-degree features, which led to its producing the model with the lowest RMSE. All linear methods converged to the same dense model. Nonlinear methods had similar performance, reaching test RMSEs between 2.0 and 2.3 rings, showing LCEN can provide interpretability without sacrificing accuracy. Finally, an LCEN model using only 3 features displayed only a 2.3% increase in RMSE relative to the best LCEN model, and another LCEN model using only 2 features displayed only a 7.1% increase in RMSE relative to the same best LCEN model. As before, LCEN built very sparse yet very accurate models. LCEN was then compared to previously published algebraic expressions and an ANN, and many other machine-learning methods, on the “Concrete Compressive Strength” dataset (Table 5). LCEN displayed a test RMSE at least 23% lower than that of the algebraic models and produces a model that is simpler and interpretable while having an RMSE competitive with that of the MLP of Ref. [34]. In the paper that originated this MLP, no type of validation is mentioned, which could indicate that the MLP results are overoptimistic. The validation performance of LCEN was better than the ANN test performance reported in Ref. [34], indicating LCEN is competitive even with ANNs on this problem if the results reported by Ref. [34] are actually validation results. Some nonlinear machine-learning architectures achieved a lower test RMSE on this task, but LCEN has the lowest validation RMSE out of all architectures tested.

Continuing these experiments, LCEN was also tested against datasets generated by human activity, as opposed to physical phenomena. In a dataset used to predict the prices of houses (“Boston housing”), LCEN achieved strong performance, with a test RMSE only 6% higher than that of an MLP, which had the best test-set performance. This exceptional result further highlights the potential of LCEN, especially in data-deficient scenarios. In the “GEFCom 2014” dynamic dataset, LCEN was able to surpass complex ANN architectures, predicting test-set power loads with a 13.1% lower RMSE than that of the former state-of-the-art model (Table 6). Due to its high robustness, LCEN also could predict for longer time horizons without large increases in the prediction error.

Overall, these experiments have demonstrated the applicability of LCEN to a multitude of scientific and nonscientific problems. LCEN models were robust to defects in the real data, including noise, multicollinearity, or sample scarcity. LCEN models were typically as accurate as or more accurate than many alternative architectures, yet were also much sparser. The LCEN architecture is also trivial to interpret and displays exactly how each input is contributing to the final output. This combination of accuracy and interpretability is essential for the deployment of machine-learning models in performance-critical scenarios, from aviation to medicine. LCEN is free, open-source, and easy to use, allowing even non-specialists in machine learning to benefit from and use it. Moreover, the additional interpretability can assist in data or model refinement efforts and can make the models robust to changes in data or adversarial input. The main limitations of LCEN are that it is not a universal function approximator, as it can model only the functions present in the expansion of dataset features, and that it sometimes is not as accurate as a dense deep learning method. If a GPU and enough time are available for model training, users in scenarios that focus on accuracy above anything else may prefer to use a deep learning architecture.

There are at least two clear future directions for this work. The first involves using the LCEN algorithm in classification tasks, as many important problems in science and engineering involve classification. A comprehensive and robust analysis of the performance of LCEN in classification tasks will follow this paper. A second future direction involves applying the LCEN algorithm to automatically generate physical equations for certain hybrid model architectures (such as physics-constrained or physics-guided ML). The potential and importance of these models has been analyzed in many works, such as in Refs. [50, 51].

The “Linear, 5-variable polynomial” dataset was created by drawing numbers from a uniform distribution between 1 and 10 and summing, potentially adding noise, such that $y=-2.8X_{0}-2.7X_{1}-5.3X_{2}+4.3X_{3}+9.0X_{4}+\epsilon$. The “Multicollinear data” dataset was created by drawing numbers from a uniform distribution between 1 and 10 to create one variable $X_{0}$, which was used together with a small amount of noise to create a correlated variable $X_{1}=X_{0}+\epsilon_{1}$; finally, they were summed such that $y=2X_{0}+2X_{1}+\epsilon_{2}$. The “Relativistic energy” dataset was created by drawing numbers from a uniform distribution between 1 and 10 or 1 and 100 for masses, and 5$\times$10⁷ and 2.5$\times$10⁸ for velocities, which represent the energy of a body as $E^{2}=c^{4}m^{2}+c^{2}m^{2}v^{2}$. With these velocity numbers, relativistic effects are responsible for 20.4% of the total squared energy on average. The “4th-degree, univariate polynomial” dataset was created by drawing numbers from a normal distribution with mean 0 and variance 5 and transforming them into the polynomial $y=X+0.5X^{2}+0.1X^{3}+0.05X^{4}+\epsilon$. The real datasets are described in Sections A6.1 and 3.2.

All possible permutations of the hyperparameters below were cross-validated.

1. 1.

   For the LASSO and Ridge regression models: $\alpha=0$ and 20 log-spaced values between -4.3 and 0 (as per np.logspace(-4.3, 0, 20)).

2. 2.

   For the elastic net (EN) models: $\alpha$ as above and L1 ratios equal to [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.97, 0.99].

3. 3.

   For the partial least squares (PLS) models: a number of components equal to all integers between 1 and a limit were used. This limit is either the number of features or 80% of the number of samples, whichever is smaller.

4. 4.

   For the LCEN models: $\alpha$ and L1 ratios as above. degree values equal to [1, 2, 3] were typically used, except when otherwise indicated (such as in the “Relativistic energy” dataset). lag = 0 was used, except for the “GEFCom 2014” dataset, which used lag = 168. cutoff values between $1$$\times$$10^{-3}$ and $5.5$$\times$$10^{-1}$ were used; higher values were used only when intentionally creating models with fewer selected features. A cutoff = 0 is used in the ablation tests for the LASSO-EN model (Section A4).

5. 5.

   For the random forest (RF) and gradient-boosted decision tree (GBDT) models: [10, 25, 50, 100, 200, 300] trees, maximum tree depth equal to [2, 3, 5, 10, 15, 20, 40], minimum fraction of samples per leaf equal to [0.01, 0.02, 0.05, 0.1], and minimum fraction of samples per tree equal to [0.1, 0.25, 0.333, 0.5, 0.667, 0.75, 1.0]. For the GBDT models, learning rates equal to [0.01, 0.05, 0.1, 0.2] were also used.

6. 6.

   For the AdaBoost (AdaB) models: [10, 25, 50, 100, 200, 300] trees/estimators and learning rates equal to [0.01, 0.05, 0.1, 0.2] were used.

7. 7.

   For the support vector machine (SVM) models: C values equal to [0.01, 0.1, 1, 10, 50, 100], epsilon values equal to [0.01, 0.025, 0.05, 0.075, 0.1, 0.15, 0.2, 0.3], and gamma values equal to [1/50, 1/10, 1/5, 1/2, 1, 2, 5, 10, 50] divided by the number of features in a dataset were used.

8. 8.

   For the multilayer perceptron (MLP), MLP with group LASSO (MLP-GL₁), and LassoNet models: the hidden layer sizes varied for each dataset. Representing an MLP with one hidden layer as [X] and an MLP with two as [X, Y], hidden layer sizes of {[800], [400], [200], [100], [800, 800], [800, 400], [400, 400], [400, 200], [200, 200], [200, 100], [100, 100]} were used with the “Diesel Freezing Point” dataset, {[18], [9], [4], [18, 18], [18, 9], [9, 9], [9, 4], [9, 2], [4, 4]} were used with the “Abalone” dataset, {[16], [8], [4], [48, 16], [48, 8], [40, 24], [40, 16], [40, 8], [32, 16], [32, 8], [24, 16], [24, 8], [16, 16], [16, 8], [8, 8], [8, 4]} were used with the “Concrete Compressive Strength” dataset, and {[26], [13], [6], [78, 26], [65, 39], [65, 26], [65, 13], [52, 39], [52, 26], [52, 13], [39, 39], [39, 26], [39, 13], [26, 26], [26, 13], [13, 13], [13, 6]} were used with the “Boston housing” dataset. Learning rates equal to [0.001, 0.005, 0.01], activation functions in [“relu”, “tanhshrink”], a batch size of 32, 100 epochs, and a cosine scheduler with a minimum learning rate equal to 1/16 of the original learning rate with 10 epochs of warm-up were also used. For the MLP-GL₁ and LassoNet, regularization parameters equal to [1$\times 10^{-4}$, 1$\times 10^{-3}$, 1$\times 10^{-2}$] were used.

To better clarify the design choices of the LCEN architecture and highlight the relevance of each individual part of the algorithm, ablation tests are performed. Three ablated architectures – LASSO-Clip (LC), EN-Clip (ENC), and LASSO-EN (LEN) – are compared with the original LCEN algorithm. Two variant architectures, LASSO-Clip-LASSO (LCL) and EN-Clip-EN (ENCEN), are also compared. The “Relativistic energy”, “Diesel Freezing Point”, “Abalone”, and “Concrete Compressive Strength” datasets are used in the ablation tests. Tests with the “Relativistic energy” dataset show that models with a Clip step had some degree of success with selecting only relevant features (Table A3). However, the ablated architectures (LC, ENC, and LEN) had much higher prediction errors for the coefficients of the relevant features, even though LC and ENC were able to select only the relevant features. The variant architectures (LCL and ENCEN) had performances closer to that of LCEN, but LCL was slightly worse in terms of error. The models that begin with EN (ENC and ENCEN) are approximately one order of magnitude slower than LCEN on all datasets (Tables A3–A6), whereas LC was approximately 50% faster than LCEN. However, LCEN consistently had the lowest validation RMSE in all datasets, and had the lowest test-set RMSE in all but one dataset. As highlighted by Table A4, LCEN built the sparsest and most accurate models out of all ablated and variant architectures trained with the “Diesel Freezing Point” dataset. Overall, these ablation experiments highlight how LCEN is the optimal architecture to maximize accuracy and selectivity while maintaining a low runtime.

All experiments were done in a personal computer equipped with a 13th Gen Intel^® Core™ i5-13600K CPU, 64 GB of DDR4 RAM, and a NVIDIA GeForce RTX 4090 GPU. Runtimes for the models in the “Diesel Freezing Point” dataset are provided in Table 4, while runtimes for LCEN and ablated architectures are provided in the tables of Section A4.