{'fraction_non_alphanumeric': 0.02998421883219358, 'fraction_numerical': 0.005322276201380078, 'mean_word_length': 4.928820399926619, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multiple Equivalent Solutions for the Lasso

Feature selection is an important problem studied in data analytics seeking to identify a minimal-size feature subset that is optimally predictive for an outcome of interest. It is also a powerful tool in Knowledge Discovery as a means for gaining domain insight, e.g., identifying which medical quantities carry unique information for the disease status. It is arguably less recognized however, that the problem may have multiple, equivalent solutions. In that case, it is misleading to domain experts to report only one of them and ignore all other equivalent solutions. In this paper, we extend a well-established single, feature selection algorithm (i.e., reporting a single solution), namely the Lasso algorithm, to the multiple solution problem based on formalized notion of equivalence for both classification and regression tasks. Empirical results are obtained using a fully automated pipeline called Just Add Data Bio or JAD Bio training and selecting multiple, linear as well as nonlinear learners, optimizing hyper-parameter values, and correcting for the bias of multiple inductions (model selection). The results show that multiple solutions do exist in real datasets, as well as the ability of the algorithm to identify a subset of them. A comparison with the Statistical Equivalent Solutions (SES) algorithm shows that Lasso equivalent solutions have better prediction performance at the cost of selecting more features.


INTRODUCTION
The problem of supervised Feature Selection (FS) (a.k.a. variable selection) has been studied for several decades in data science fields, such as statistics, machine learning, and data mining. Informally, the problem can be defined as selecting a subset of the available feature set, such that it is of minimal-size and at the same time maximally-predictive for an outcome of interest. Feature selection is performed for several reasons. It may reduce the cost or risk of measuring, computing, storing, and processing of the features. It leads to models of smaller dimensionality that are easier to visualize, inspect, and understand. It may even result in more accurate models by removing the noise, treating the curse-of-dimensionality, and facilitating model fitting. Arguably however, feature selection is primarily employed as a means for knowledge discovery and gaining intuition into the mechanisms generating the data. Indeed, a deep theoretical connection with causality has been identified [1]. Actually, it is often the case that FS is the primary goal of an analysis and the learned model is only a by-product. For example, a medical researcher analyzing their molecular data may not care about the classification model to cancerous or healthy tissue; they can diagnose the tissues themselves without the use of the computer on the microscope. The whole point of their analysis is to identify the features (e.g., gene expressions) that carry all the information to perform the diagnosis.
When FS is employed for knowledge discovery and understanding considering multiple, equivalent solutions is paramount. It is misleading to present to the domain expert a single feature subset and proclaim all other features are either redundant or irrelevant, if there are multiple equivalent solutions. In addition, when features have an associated measurement cost, one seeks the least cost solution. We argue that in these cases, the analyst ought to identify all solutions and present the possibilities to the domain expert. Multiple, equivalent FS solutions often exist in practice [2], especially in fields with low sample size, high dimensionality and noisy features. Biology and medicine are prototypical fields where these adverse conditions are met. In a seminal paper [3], Ein-dor and co-authors demonstrate that multiple, equivalent prognostic signature for breast cancer can be found just by analyzing the same dataset with a different partition in training and test set, showing that several genes exist which are practically exchangeable in terms of predictive power. Statnikov and Aliferis [4] further show that the presence of multiple optimal signatures is not a rare occurrence, and is actually common in biological datasets.
Despite the relevance of the problem, to date only few algorithms address it directly. Two constraint-based approaches (constraints stemming from the results of conditionally independence tests) are the Target Information Equivalence (TIE * [5]) and the Statistically Equivalent Signatures (SES [2]) algorithms. Another algorithm which has been recently proposed by Cox and Battey [6] searches for multiple solutions across a large number of separate analyses. Unfortunately, both the Cox-Battey method and TIE * are quite computationally intensive.
In this work, we show that multiple equivalent solutions with controlled performance deviations can be also defined and identified for the popular least absolute shrinkage and selection operator (Lasso) method [7]. Previous studies based on bootstrapping [8], sub-sampling [9] or locally perturbing the support of the Lasso solution [10] propose to run the Lasso multiple times and collectively establish the best (unique) solution. Their goal is to improve the ordinary Lasso solution in terms of either accuracy or stability. However, this is irrelevant to our goal which is to determine the set of equivalent Lasso solutions.
The case when strongly-equivalent Lasso solutions (SELSs) exist has been developed and studied in [11]. The author provided there various characterizations of the space of all SELSs. Multiple SELSs however, rarely occur in real datasets. Henceforth, we propose a relaxed definition of equivalence requiring the equivalent solutions to have not exactly the same, but a similar (within a tolerance threshold) in-sample loss. This is useful because with finite sample different feature subsets may appear statistically indistinguishable in terms of loss, even though asymptotically there is a single optimal solution. Specifically, using the root mean squared error (RMSE) as the loss function, we define the RMSE-equivalent Lasso solutions (RELSs). In contrast, a natural loss for l 1 -penalized logistic regression (i.e., logistic Lasso) is the deviance leading to Deviance-equivalent Lasso solutions (DELSs). Inspired by the characterization of the SELS space, we devise an algorithm that computes a maximal subspace of the space of relaxed Lasso solutions. A set of constraints that guarantees sparsity as well as the key characteristics of the Lasso solution(s) is imposed. The proposed algorithm handles both RMSE and deviance equivalence in a unified manner and actually it is easily extensible to any convex performance metric. Furthermore, we derive theoretical bounds both for the RMSE and the deviance that relate the allowed tolerance with the spectral properties of the predictors matrix.
We perform empirical experiments on several real datasets providing corroborating evidence of the prevalence of multiple, equivalent solutions in real scenarios and quantifying the efficacy of the proposed Lasso extension. Indeed, the presence of multiple heterogeneous signatures (i.e., solutions) having at the same time low value for their coefficient of variation indicates that multiple, equivalent solutions are common across several application fields. In addition, we compare the proposed algorithm against the SES algorithm that can scale w.r.t. the dimensionality of the data. A prominent difference is that Lasso equivalent solutions have on average better prediction performance at the cost of selecting larger signatures while the SES algorithm is more parsimonious. A novelty of the comparison is the employment of a fully automated analysis pipeline (called Just Add Data Bio or JAD Bio) that tries both linear and non-linear learners, optimizes their hyper-parameter values, and corrects the estimation of performance for the bias of the tuning. Thus, the results are obtained after matching each feature selection with the most appropriate learner and its hyper-parameter values.

DEFINITION OF LASSO SOLUTIONS EQUIVALENCE 2.1 Lasso preliminaries
We first introduce the Lasso inference problem. Given an outcome vector (a.k.a. target variable) y ∈ R n , a matrix X ∈ R n×p with the predictor variables and the coefficient vector β ∈ R p , the Lagrangian form of Lasso is defined as where λ ≥ 0 is a penalty parameter. Lasso solvers such as LARS [12] and FISTA [13] have been extensively applied in thousands of real problems while extensions of (1) to generalized linear models (GLMs) such as the logistic regression with l 1 -norm regularization have been also proposed [14]. In GLM Lasso, the quadratic cost term is replaced by the deviance which essentially behaves as the negative log-likelihood.

Strong equivalence
We briefly present the strong equivalence and its key properties which have been firstly defined and developed in [11]. It has been proved that when the entries of X are drawn from a continuous probability distribution, the uniqueness of the Lasso solution is almost surely guaranteed [11]. However, when the predictor variables take discrete values or there exists perfect collinearity then several optimal Lasso solutions may exist. Mathematically, given λ > 0, two vectorsβ,β ′ ∈ R p are Lasso equivalent solutions in the strong sense if and only if L(β) = L(β ′ ). The set of all SELS is defined as Additionally, two SELSs not only share the same cost function but also predict the same values for the target variable [11] (i.e., Xβ = Xβ ′ ). Consequently, SELSs have the same l 1 norm, ||β || 1 = ||β ′ || 1 .
Another characteristic of all SELSs is that the non-zero coefficients are not allowed to flip their sign [11]. Hence, ifβ i > 0 for the i-th coefficient of a Lasso solution thenβ ′ i ≥ 0 for any SELS,β ′ , and, similarly, for the coefficients with negative values.
Proceeding, letβ be the Lasso solution with the largest active set (or support), E := supp(β) := {i :β i 0} which is also known as the equicorrelation set [11]. The maximal Lasso solution is computed using either the elastic net with vanishing l 2 penalty norm or the variant of LARS described in [11]. With a slight abuse of notation due to the restriction of the solution space to the predictor variables with non-zero coefficients, the set of all SELSs can be rewritten as where X E is the matrix that contains only the columns of X that are indexed by the equicorrelation set E while S := diaд(s) is the diagonal matrix with the signs of the Lasso solution, s := sign(β E ).
The first constraint ensures that all elements in K will have the same fitted value while the second constraint ensures that coefficients' sign will not flip.

Enumeration Algorithm for SELSs.
A bounded polyhedron (i.e., a polytope) can be represented either as the convex hull of a finite set of vertices or by using a combination of linear constraint equalities and inequalities. In particular, the vertices of the convex polytope, K, defined above corresponds to the "extreme" SELSs. Thus, enumerating the vertices of K from the set of equality and inequality constraints, we can enumerate all SELSs. There exist algorithms that enumerate the vertices defined by a set of inequality constraints (see Fukuda et al. [15] and the references therein) and they can be extended to take into account equality constraints, too. In this paper, we employed the Matlab package by Jacobson [16] which contains tools for converting between the (in)equality and the vertices representations.

Cardinality of equivalent solutions.
It is noteworthy that the number of vertices, like the number of edges and faces, can grow exponentially fast with the dimension of the polytope making the enumeration algorithm impractical when the dimension is larger than 20. In such cases, we cannot enumerate all the extreme SELSs in reasonable time. Nevertheless, it would be useful to know which of the variables participate in all SELSs and which are not. In Section 4, a practical categorization of the variables into dispensable (participate in some solutions) and indispensable (participate in all solutions) which has been firstly introduced in [11] is presented and extended to the relaxed equivalence, too.

Relaxed equivalence
A less restrictive equivalence is to require two solutions to have similar performance (loss, error). We will search for solutions whose performance metric differ by a small tolerance from a given solution. Moreover, in order to avoid handling absolute quantities, we suggest working with the relative performance. Denoting by D(·) the performance metric,β the given solution and TOL the tolerance, we say thatβ is performance equivalent toβ if and only if the following relation on the relative performance metric is satisfied However, controlling only the performance metric may add unnecessary redundancy to the relaxed solutions destroying the desired property of sparsity. Indeed, a potentially large number of irrelevant or redundant variables may satisfy the performance constraint for a given TOL, nonetheless, they should not belong to the set of equivalent solutions. To alleviate this issue, we propose to restrict the space of relaxed solutions by allowing only the non-zero coefficients of the given solution to vary. Thus, for a given Lasso solution,β with support supp(β) = {i :β i 0}, a performance metric, D(·), and a tolerance TOL, we define the set of D-equivalent Lasso solutions as where −A denotes the complement set of A. The following proposition states that the convexity of the performance metric with respect to β is inherited to the set of all relaxed Lasso solutions. Proposition 1. Let D(·) be a convex performance metric. Then, and convexity is proved. ■ Relaxed equivalence can be defined for several different performance metrics, which depend on the type of the analysis task (classification, regression, survival analysis, etc.). Different characteristics of the relaxed solution space can be conveyed by a suitable performance metric. For instance, when the Lasso cost function is used as a performance metric, the support ofβ is maximal and TOL is set to 0 then the strong equivalence is obtained. Next, we present one performance metric for regression problems and one for classification tasks.
is a standard performance metric for regression problems. The set of RMSE-equivalent Lasso Solutions (RELSs), K T O L RMS E (β), is a convex set since the RMSE is a convex performance metric. We note also that RMSE-equivalence with tolerance TOL is almost equivalent with MSE-equivalence with tolerance 2TOL since RMSE(β) for small values of the tolerance.
2.3.2 Deviance equivalence for logistic Lasso. Logistic Lasso is preferred as a feature selection method in binary classification tasks compared to the ordinary Lasso because it takes into consideration the distinguishing properties of the problem. In logistic Lasso, the quadratic term in (1) is replaced by the deviance. Hence, despite the fact that RMSE is applicable, a more appropriate performance metric is the mean deviance defined as The set of Deviance-equivalent Lasso Solutions (DELSs), K T O L DEV (β), is also a convex set since the deviance is a convex function of its argument.

ENUMERATING RELAXED LASSO SOLUTIONS
In the definition of relaxed equivalence (i.e., (5)), we limit the active set of relaxed equivalent solutions to a subset of the given solution's active set. Despite this restriction, the enumeration of all solutions is still difficult because of the shape of K T O L D (β); the set continuously curves resulting in an infinite number of extreme solutions. Hence, we propose to restrict the solutions that we will eventually enumerate to a subset of K T O L D (β). Inspired by the representations of SELSs, we propose an approach that enumerates the vertices of the largest convex polytope of the relaxed space that allows relative performance be less than TOL. The central idea is to relax the Null space constraint while keeping intact the constraint on the sign of the coefficients. For this purpose, we rewrite the Null space constraint using the singular value decomposition for X E . Let X E be decomposed as where U E and V E are orthogonal matrices with the left and right space eigenvectors, respectively, while Σ E is a diagonal matrix with diagonal elements the ordered singular values (i.e., corresponds to the matrix with the eigenvectors whose singular values are non-zero.
Inspired by the above representation, we propose to relax the constraint to as above while i * is an integer between 1 and |E | to be specified later. Thus, the convex polytope for the relaxed Lasso solutions is defined as We further constrain the relaxed Lasso solutions to live in a box and have coefficients that are at most l far away from the given Lasso solution,β. This constraint is added so as to guarantee the boundedness of the polyhedron defined by the other two constraint. We choose to set the box size to be l = ||β || ∞ = max iβi which allows any coefficient to vary till the zero value. Since the relaxed polytope is performance metric ignorant, our goal is to choose appropriately i * so as K * is a maximal subset of the relaxed solution space, K T O L D (β). Figure 1 demonstrates a two-dimensional example where the Lasso regression problem has a unique solution (small red circle). By setting i * = 1 the relaxed set of Lasso solutions constitutes the red line that connects the positive axes determined by the direction of v 2 . The two vertices (solid red dots) are candidates as "extreme" relaxed Lasso solutions which are further tested for ensuring that their relative performance is below the maximum tolerance. Before proceeding with an algorithm that enumerates the vertices of K * that is guaranteed to be a subset of the relaxed space, we derive theoretical bounds for two performance metrics.

Theoretical Bounds
(b) Letβ be a logistic Lasso solution. For anyβ ∈ K * it holds that

Enumeration algorithm
Algorithm 1 reliably enumerates relaxed solutions that belong to a subset of the set of all relaxed Lasso solutions. It takes as input a Lasso solution, a convex performance metric D(·), a tolerance TOL and the maximum dimension d max of the subspace that will be explored. We assume that the provided Lasso solution has maximal support hence the equicorrelation set and the sign vector are directly estimated as in lines 2 and 3 of Algorithm 1. In the for loop, the dimension of the subspace is increased starting from 1 until the maximum value is reached or there are vertices whose performance metric is above the tolerance. In other words, the algorithm removes the restriction stemming from the component v with the least singular value, allowing itself to search along the v direction. It sequentially removes restrictions corresponding to components with larger singular values. As the dimension of the exploration subspace is increased the possibility of finding solutions whose performance exceeds the tolerance is also increased. If such solutions occur then the algorithm breaks and the solutions that exceed the tolerance are discarded. Due to the convexity of the performance metric, the polytope that is defined by the remaining solutions is also a subset of the relaxed space. Thus, the algorithm is sound. Additionally, a variant of binary search can be utilized instead of the linear search for identifying the optimal i * . However, the performance gains in computational time is minimal due to the fact that the number of vertices grows exponentially fast since the computational cost for the last iteration is proportional to the total cost of all previous iterations.  return only the D-equivalentβ's Remark. We also investigated alternative definitions of the subset, K * , where we relax different properties of the stronglyequivalent characterization. For instance, we relaxed the condition X (β −β) = 0 to ||X (β −β)|| <= r with r > 0 resulting, however, in optimization problems which are algorithmically intractable. Overall, the choice of performance measures (RMSE and Deviance) and the particular subset, K * , in (7), were based on the mathematical simplicity and the practical implementation both being very critical for the adoption of a less popular but important idea of searching for and enumerating multiple equivalent solutions.

Determining the reference Lasso solution
Algorithm 1 searches for equivalent solutions in the least sensitive directions of the performance metric. It may discard features from the reference support but it cannot include features not initially selected. Therefore, it is important to provide an initial solution with the maximum support. We compute the reference Lasso solutionβ, by optimally tuning its penalty parameter. We perform cross-validated hyper-parameter tuning using AUC (resp. MAE) as performance metric for classification (resp. regression) problems. It has been highlighted and proved [17, Proposition 1] that many noise features are included in the prediction-oracle solution (i.e., the reference solution we estimate). Indeed, the value of the penalty parameter is lower in the prediction-oracle solution than the consistent one hinting towards larger support. We also observe similar behavior in our experiments with real datasets in the Results section. Thus, we expect the respective initial support to contain most, if not all, of the related to the target variable features. Of course, alternative methods based on elastic net or bootstrapping for selecting the reference solution/support can be utilized, however, such an exploration is beyond this paper's scope and it is left as future work.

VARIABLE CATEGORIZATION
Due to the exponential growth of the polytope's vertices, the complete enumeration is not always feasible. Nevertheless, there is an alternative approach to qualitatively assess the predictor variables by summarizing the interesting findings, e.g., reporting the set of features that belong to all solutions or only to some. Exploiting the fact that the feasible range of a coefficient's value for all Lasso solutions can be efficiently computed as it was shown in [11], a variable categorization is possible. We present the existing formulation for SELSs and based on it we extend it for RELSs. SELSs. For each i ∈ E, the i-th coefficient's lower boundβ l i and upper boundβ u i are computed by solving the linear programŝ respectively. If 0 is an element of the set [β l i ,β u i ] then the i-th variable is called dispensable otherwise it is called indispensable and they participate in all Lasso solutions. The fact that the sign of a coefficient remains the same in all SELSs implies that dispensable variables have eitherβ l = 0 orβ u = 0 while indispensable variables have eitherβ l > 0 orβ u < 0. From a practical perspective, the number of linear programs to be solved is 2|E | which is feasible. Relaxed equivalence. Similarly, we can relax the criterion on dispensable/indispensable variables. Since the computational cost of the linear programs that are solved is manageable, we can discard the constraint on the control of the subspace dimension (i.e., there is no need for d max ). The linear programs for the lower and upper bounds for a given i * arē respectively. Evidently, as we increase i * the number of dispensable variables is increased while the number of indispensable variables is decreased. While in this paper we follow the terminology of "dispensable" and "indispensable" variables from [11], we would like to emphasize that it may be misleading. It implies that dispensable variables could be dispensed with and filtered out in some sense. Consider however, two featuresW and Z that are perfect copies of each other and at the same time, the strongest predictors. According to the definition they are dispensable but neither can be filtered out without a measurable loss in performance. Instead, we have recently proposed the terms "indispensable" and "replaceable" 2 to clearly indicate that they the latter features can be replaced, not filtered out.

Just Add Data Bio Engine
JAD Bio employs a fully-automated machine learning pipeline for producing a predictive model given a training dataset, and an estimate of its predictive performance. JAD Bio performs multiple feature selection and predictive modeling. It tries several configurations, i.e., combinations of preprocessing algorithms, feature selection algorithm, predictive modeling algorithms, and values of their hyper-parameters using a grid-search in the space of hyperparameters. The best configuration is determined using K-fold Cross Validation, which is then used to produce the final models for each feature subset on all data. The cross-validated performance of the best configuration is known to be optimistic due to the multiple tries [18]; JAD Bio estimates and removes the optimism using a bootstrap method before returning the final performance estimate [19]. The final performance estimates of JAD Bio are actually slightly conservative. Specifically in terms of algorithms, JAD Bio trains several basic and advanced, linear and non-linear, multivariate machine learning and statistical models; namely, for classification problems trains Support Vector Machine models (SVMs) with linear, full polynomial, and Gaussian kernels, Ridge Logistic Regression models, Random Forests models, and Decision Trees, while for regression problems trains Ridge Linear Regression models instead of the logistic ones. In the following set of experiments more than 500 configurations were tried each time to find the optimal ones, resulting in more than 10000 models trained to compare two feature selection algorithms over nine datasets. JAD Bio is able to incorporate other user-defined preprocessing, transformation, and feature selection algorithms written in Java, R, or Matlab programming languages. This functionality allowed us to incorporate the proposed feature selection method into the pipeline.

Evaluation Pipeline and Tuning
Two variants of JAD Bio were used and compared: the original with the SES algorithm for the feature selection and an alternative with a Matlab implementation of the proposed method. To evaluate and compare them we split each dataset in two stratified parts. Then, JAD Bio is applied on each half (as training set) for obtaining a set of predictive models and their performance estimation, corresponding to different feature subsets (i.e., solutions found by the feature selection algorithms). The trained models are then tested on the remaining data (as test set) to compute the performance on new data. We then test whether the performance of the solutions found are actually similar. Since we partition the data into two, equivalence of the performance of the feature selection subset is determined using the same sample size as in training and thus having the same statistical power. The new proposed method was tuned with penalty parameter, λ, taking values from 10 −3 to 10 3 with multiplicative step 10 0.1 . We set the TOL value to 0.01 allowing Lasso solutions having up to 1% performance error in the training set.

RESULTS
6.1 Multiple, equivalent solutions are common for predictive analytics tasks Table 2 reports the Lasso and SES results. Both algorithms identify multiple solutions on almost all datasets and splits, indicating that the presence of multiple solutions is a common problem across several application domains. At the same time, the coefficient of variation (CoV, defined as the ratio between standard deviation and average) for the hold-out performances is always quite low; the median CoV for Lasso is 0.0097 and for SES it is 0.0049 (maximum values are respectively 0.0421 and 0.0373). This means that the different solutions produce models with performance having ∼ 0.01 standard deviation around the mean AUC value. These results support the claim that multiple solutions found by Lasso and SES are indeed equivalent in terms of predictive power.
Interestingly, when SES retrieves a large number of signatures (i.e., solutions) so does Lasso, and vice versa. Both algorithms provide numerous solutions for the arcene and secom datasets, only a few solutions for Ovarian, Parkinson, prostate, vV, BC_Continuous, and MP (an exception is the bankruptcy dataset). This is despite the fact that the two algorithms follow quite different approaches. Thus, the results provide evidence that the order of magnitude for the number of equivalent solutions is a characteristic of the data at hand. Figure 2 graphically represents the hold-out performances (x-axis, AUC for classification datasets and R 2 for regression ones) and number of selected features (y-axis) for each signature retrieved by Lasso (circles) and SES (triangles) for the first split of each dataset (distinguished by color). The figure shows that Lasso consistently selects up to one order of magnitude more variables than SES (Table 2). In nine out of 18 comparisons (9 datasets × 2 splits each) Lasso also achieves better performances than SES (adjusted p-value ≤ 0.05), with an average AUC difference of 0.034 between the Table 2: Lasso and SES comparison. For each dataset, split and method (Lasso/SES), the table reports the training set sample size, the number of identified signatures, the coefficient of variation (CoV) across the multiple signatures for the hold-out performance (AUC for classification, R 2 for regression), the CoV for the number of selected variables, the average hold-out performance (along with the two-tail, t-test adjusted p-value assessing the difference between SES and Lasso), as well as the average number of selected variables (along with the two-tail, t-test adjusted p-value assessing that the average number of variable selected by Lasso is different from the fixed number of variables selected by SES). Classification datasets are reported first, while the regression ones are at the bottom.  The Jaccard index between two sets is defined as the ratio between the size of their intersection and the size of their union, and ranges from 0 (disjoint sets) to 1 (perfectly overlapping sets), with two equal-size sets sharing half of their elements achieving 0.33. For computing the number of solutionspecific features between two sets we first take their union and then subtract their intersection. The resulting features are either in one solution or in the other, but not in both. For the cases where numerous solutions are selected by SES (arcene and secom datasets), both the Jaccard index and the number of dispensable elements indicate a quite elevated heterogeneity, meaning that the multiple solutions proposed by SES include solutions quite different from each other. When only a handful of signatures are retrieved by SES, they tend to be more homogeneous. This is a consequence of the mechanisms used by SES for producing multiple signatures, which constraints each signature to differ from its closest sibling only for one variable. In contrast, Lasso solutions can have more than two signature-specific features even when only two signatures are retrieved. This can be seen in several rows of Table 3, e.g., for the arcene dataset (signature size 53, 54) and for the secom dataset (signature size 43). Interestingly, in several cases Lasso solutions differ from each other on average by 5 or more solutionspecific features, confirming that the multiple solutions retrieved by Lasso differ from each other due to the joint replacement of several variables.

CONCLUSION
In this paper we present an algorithm for deriving multiple, equivalent Lasso solutions with theoretical guarantees. The algorithm is empirically evaluated using an automatic platform, JAD Bio, over several learners and tuning hyper-parameters; the same platform is employed for a comparative evaluation against the SES algorithm. The results confirm that (a) multiple Lasso solutions are often found in real-world datasets, (b) these multiple solutions can substantially differ from each other, despite their similar performances, and (c) Lasso and SES provide comparable results, with the first being on average slightly more predictive at the cost of selecting a larger number of features.