Utilizing evolutionary algorithms for quick function choice with massive datasets

11 min learn

·

12 hours in the past

That is half 1 of a two-part sequence about function choice. Learn half 2 right here.

If you’re becoming a mannequin to a dataset, you might have to carry out function choice: conserving just some subset of the options to suit the mannequin, whereas discarding the remainder. This may be vital for a wide range of causes:

• to maintain the mannequin explainable (having too many options makes interpretation more durable)
• to keep away from the curse of dimensionality
• to maximise/reduce some goal perform associated to the mannequin (R-squared, AIC, and many others)
• to keep away from a poor match, and many others.

If the variety of options N is small, an exhaustive search could also be doable — you possibly can actually attempt all potential mixtures of options, and simply maintain the one which minimizes the associated fee / goal perform. But when N is massive, then an exhaustive search may be inconceivable. The whole variety of mixtures to attempt is 2^N which, if N is bigger than just a few dozen, turns into prohibitive (it’s an exponential perform). In such instances you should use a heuristic technique: discover the search area in an environment friendly means, in search of a mix of options that minimizes the target perform you utilize to conduct the search.

What you’re in search of is a vector [1, 0, 0, 1, 0, 1, 1, 1, 0, ...] of size N, with parts taking values from {0, 1} . Every ingredient within the vector is assigned to a function; if the ingredient is 1, the function is chosen; if the ingredient is 0, the function is discarded. You’ll want to discover the vector that minimizes the target perform. The search area has as many dimensions N as there are options, and the one potential values alongside any dimension are 0 and 1.

Discovering a superb heuristic algorithm is just not trivial. The regsubsets() perform in R has some choices you should use. Additionally, scikit-learn provides a number of strategies to carry out a heuristic function choice, offered your drawback is an efficient match for his or her strategies. However discovering a superb, basic objective heuristic — in essentially the most basic kind — is a tough drawback. On this sequence of articles we’ll discover just a few choices that will work moderately effectively even when N is sort of massive, and the target perform may be actually something, so long as it’s straightforward to compute.

For all optimization strategies on this sequence of articles, I’m utilizing the highly regarded Home Costs dataset on Kaggle (MIT license) — after just a few easy function transformations, it finally ends up having 213 options (N=213) and 1453 observations. The mannequin I’m utilizing is linear regression, statsmodels.api.OLS() , and the target perform I’m attempting to attenuate is BIC — the Bayesian Info Criterion, a measure of knowledge loss, so a decrease BIC is healthier. It’s much like AIC — the Akaike Info Criterion — besides BIC tends to provide extra parsimonious fashions (it prefers fashions with fewer options). Minimizing both AIC or BIC tends to scale back overfitting. However different goal capabilities may be used as an alternative, e.g. R-squared (the defined variance within the goal) or adjusted R-squared — simply remember that a bigger R-squared is healthier, in order that’s a maximization drawback.

Finally, the selection of goal perform is irrelevant right here. What issues is that we have now one goal perform that we’re constantly attempting to optimize utilizing varied strategies.

The complete code used on this sequence of articles is contained in a single pocket book in my function choice repository — additionally linked on the finish. I’ll present code snippets throughout the textual content right here, however please examine the pocket book for the total context.

However first, let’s examine a well known, tried-and-true function choice approach, which we’ll use as a baseline comparability with the extra advanced strategies described later.

SFS, the ahead model, is a straightforward algorithm. It begins by attempting to decide on a single function, and it selects that function that minimizes the target perform. As soon as a function is chosen, it stays chosen eternally. Then it tries so as to add one other function to it (for a complete of two options), in such a means as to attenuate the target once more. It will increase the variety of chosen options by 1, each time looking for one of the best new function so as to add to the present assortment. The search ends when all options are tried collectively. Whichever mixture minimizes the target, wins.

SFS is a grasping algorithm — every alternative is domestically optimum — and it by no means goes again to right its personal errors. But it surely’s moderately quick even when N is sort of massive. The whole variety of mixtures it tries is N(N+1)/2 which is a quadratic polynomial (whereas an exhaustive search must carry out an exponential variety of trials).

Let’s see what the SFS code might appear to be in Python, utilizing the mlxtend library:

import statsmodels.api as sm
from mlxtend.feature_selection import SequentialFeatureSelector as SFS
from sklearn.base import BaseEstimator
class DummyEstimator(BaseEstimator):
# mlxtend desires to make use of an sklearn estimator, which isn't wanted right here
# (statsmodels OLS is used as an alternative)
# create a dummy estimator to pacify mlxtend
def match(self, X, y=None, **kwargs):
return self
def neg_bic(m, X, y):
# goal perform
lin_mod_res = sm.OLS(y, X, hasconst=True).match()
return -lin_mod_res.bic
seq_selector = SFS(
DummyEstimator(),
k_features=(1, X.form[1]),
ahead=True,
floating=False,
scoring=neg_bic,
cv=None,
n_jobs=-1,
verbose=0,
# make sure that the intercept is just not dropped
fixed_features=['const'],
)
n_features = X.form[1] - 1
objective_runs_sfs = spherical(n_features * (n_features + 1) / 2)
t_start_seq = time.time()
# mlxtend will mess along with your dataframes in case you do not .copy()
seq_res = seq_selector.match(X.copy(), y.copy())
t_end_seq = time.time()
run_time_seq = t_end_seq - t_start_seq
seq_metrics = seq_res.get_metric_dict()

It runs by way of the mixtures shortly and these are the abstract outcomes:

greatest okay:         36
greatest goal: 33708.98602877906
R2 @ greatest okay:    0.9075677543649224
goal runs: 22791
complete run time: 44.850 sec

The very best variety of options is 36 out of 213. The very best BIC is 33708.986, and it takes lower than 1 minute to finish on my system. It invokes the target perform 22.8k instances.

These are one of the best BIC and R-squared values, as capabilities of the variety of options chosen:

For extra data, such because the names of the options truly chosen, examine the pocket book within the repository.

Now let’s attempt one thing extra advanced.

This can be a numeric optimization algorithm. It’s in the identical class as genetic algorithms (they’re each evolutionary), however CMA-ES is sort of distinct from GA. It’s a stochastic algorithm, and it’s derivative-free — it doesn’t require you to compute a by-product of the target perform (not like gradient descent, which depends on partial derivatives). It’s computationally environment friendly, and it’s utilized in a wide range of numeric optimization libraries corresponding to Optuna. I’ll solely try a short sketch of CMA-ES right here; for extra detailed explanations, please seek advice from the literature within the hyperlinks part on the finish.

Contemplate the 2-dimensional Rastrigin perform:

The heatmap beneath exhibits the values of this perform — brighter colours imply the next worth. The perform has the worldwide minimal within the origin (0, 0), however it’s peppered with many native extremes. We wish to discover the worldwide minimal by way of CMA-ES.

CMA-ES is predicated on the multivariate regular distribution. It generates take a look at factors within the search area from this distribution. You’ll have to guess the unique imply worth of the distribution, and its commonplace deviation, however after that the algorithm will iteratively modify all these parameters, transferring the distribution within the search area, in search of one of the best goal perform values. Right here’s the unique distribution that the take a look at factors are drawn from:

xi is the set of factors that the algorithm generates at every step, within the search area. lambda is the variety of factors generated. The imply worth of the distribution will probably be up to date at each step, and hopefully will converge finally on the true answer. sigma is the usual deviation of the distribution — the unfold of the take a look at factors. C is the covariance matrix: it defines the form of the distribution. Relying on the values of C the distribution might have a “spherical” form or a extra elongated, oval form. Modifications to C enable CMA-ES to “sneak” into sure areas within the search area, or keep away from different areas.

A inhabitants of 6 factors was generated within the picture above, which is the default inhabitants dimension picked by the optimizer for this drawback. This is step one. After this, the algorithm must:

• compute the target perform (Rastrigin) for every level
• replace the imply, the usual deviation, and the covariance matrix, successfully creating a brand new multivariate regular distribution based mostly on what it has realized from the target perform
• generate a brand new set of take a look at factors from the brand new distribution
• repeat till some criterion is fulfilled (both converge on some imply worth, or exceed the utmost variety of steps, and many others.)

I cannot present the updates for all distribution parameters right here, or else this text will turn out to be very massive — examine the hyperlinks on the finish for an entire clarification. However updating simply the imply of the distribution is sort of easy, and it really works like this: after computing the target perform for every take a look at level, weights are assigned to the factors, with the bigger weights given to the factors with a greater goal worth, and a weighted sum is computed from their positions, which turns into the brand new imply. Successfully, CMA-ES strikes the distribution imply in direction of the factors with a greater goal worth:

If the algorithm converges to the true answer, then the imply worth of the distribution will converge to that answer. The usual deviation will converge to 0. The covariance matrix will change the form of the distribution (spherical or oval), relying on the geography of the target perform, extending into promising areas, and shying away from dangerous areas.

Right here’s an animated GIF displaying the evolution in time of the take a look at factors whereas CMA-ES is fixing the Rastrigin drawback:

The 2D Rastrigin perform was comparatively easy, because it solely has 2 dimensions. For our function choice drawback, we have now N=213 dimensions. Furthermore, the area is just not steady. Every take a look at level is an N-dimensional vector with element values from {0, 1} . In different phrases, every take a look at level appears like this: [1, 0, 0, 1, 1, 1, 0, ...] — a binary vector. However aside from that, the issue is identical: we have to discover the factors (the vectors) that reduce the target perform: the BIC parameter of an OLS mannequin.

Right here’s a easy model of the CMA-ES code for function choice, utilizing the cmaes library:

def cma_objective(fs):
features_use = ['const'] + [
f for i, f in enumerate(features_select) if fs[i,] == 1
]
lin_mod = sm.OLS(y_cmaes, X_cmaes[features_use], hasconst=True).match()
return lin_mod.bic
X_cmaes = X.copy()
y_cmaes = y.copy()
features_select = [f for f in X_cmaes.columns if f != 'const']
dim = len(features_select)
bounds = np.tile([0, 1], (dim, 1))
steps = np.ones(dim)
optimizer = CMAwM(
imply=np.full(dim, 0.5),
sigma=1 / 6,
bounds=bounds,
steps=steps,
n_max_resampling=10 * dim,
seed=0,
)
max_gen = 100
best_objective_cmaes_small = np.inf
best_sol_raw_cmaes_small = None
for gen in tqdm(vary(max_gen)):
options = []
for _ in vary(optimizer.population_size):
x_for_eval, x_for_tell = optimizer.ask()
worth = cma_objective(x_for_eval)
options.append((x_for_tell, worth))
if worth < best_objective_cmaes_small:
best_objective_cmaes_small = worth
best_sol_raw_cmaes_small = x_for_eval
optimizer.inform(options)
best_features_cmaes_small = [
features_select[i]
for i, val in enumerate(best_sol_raw_cmaes_small.tolist())
if val == 1.0
]
print(f'greatest goal: {best_objective_cmaes_small}')
print(f'greatest options:  {best_features_cmaes_small}')

The CMA-ES optimizer is outlined with some preliminary guesses for the imply worth and for sigma (commonplace deviation). It then loops by way of many generations, creating take a look at factors x_for_eval , evaluating them with the target, modifying the distribution (imply, sigma, covariance matrix), and many others. Every x_for_eval level is a binary vector [1, 1, 1, 0, 0, 1, ...] used to pick out the options from the dataset.

Please notice that the CMAwM() optimizer is used (CMA with margin) as an alternative of the default CMA(). The default optimizer works effectively for normal, steady issues, however the search area right here is high-dimensional, and solely two discrete values (0 and 1) are allowed. The default optimizer will get caught on this area. CMAwM() enlarges a bit the search area (whereas the options it returns are nonetheless binary vectors), and that appears sufficient to unblock it.

This straightforward code positively works, however it’s removed from optimum. Within the companion pocket book, I’ve a extra advanced, optimized model of it, which is ready to discover higher options, quicker. However the code is huge, so I cannot present it right here — examine the pocket book.

The picture beneath exhibits the historical past of the advanced, optimized CMA-ES code trying to find one of the best answer. The heatmap exhibits the prevalence / reputation of every function at every technology (brighter = extra widespread). You may see how some options are all the time highly regarded, whereas others fall out of style shortly, and but others are “found” later. The inhabitants dimension picked by the optimizer, given the parameters of this drawback, is 20 factors (people), so function reputation is averaged throughout these 20 factors.