Utilizzo di algoritmi evolutivi per la selezione rapida delle funzionalità con set di dati di grandi dimensioni
Questa è la prima parte di una serie in due parti sulla selezione delle funzionalità. Leggere parte 2 qui.
Quando adatti un modello a un set di dati, potrebbe essere necessario eseguire la selezione delle funzionalità: mantenendo solo alcuni sottoinsiemi di funzionalità per adattarli al modello, scartando il resto. Ciò può essere necessario per diversi motivi:
- mantenere il modello spiegabile (avere troppe caratteristiche rende l’interpretazione più difficile)
- per evitare la maledizione della dimensionalità
- massimizzare/minimizzare alcune funzioni obiettivo legate al modello (R quadrato, AIC, ecc.)
- per evitare un adattamento inadeguato, ecc.
Se il numero di funzionalità N è piccolo, può essere fattibile una ricerca esaustiva: puoi letteralmente provare tutte le possibili combinazioni di funzionalità e mantenere solo quella che minimizza la funzione costo/obiettivo. Ma se N è grande, una ricerca esaustiva potrebbe essere impossibile. Il numero totale di combinazioni da provare è 2^N la quale, se N è maggiore di qualche decina, diventa proibitiva (è una funzione esponenziale). In questi casi è necessario utilizzare un metodo euristico: esplorare lo spazio di ricerca in modo efficiente, cercando una combinazione di caratteristiche che riduca al minimo la funzione obiettivo utilizzata per condurre la ricerca.
Quello che stai cercando è un vettore (1, 0, 0, 1, 0, 1, 1, 1, 0, ...) di lunghezza N, con elementi che prendono valori da {0, 1} . Ogni elemento nel vettore è assegnato a una caratteristica; se l’elemento è 1 la caratteristica è selezionata; se l’elemento è 0, la funzionalità viene scartata. Devi trovare il vettore che minimizza la funzione obiettivo. Lo spazio di ricerca ha tante dimensioni N quante sono le caratteristiche e gli unici valori possibili lungo qualsiasi dimensione sono 0 e 1.
Trovare un buon algoritmo euristico non è banale. IL regsubsets() la funzione in R ha alcune opzioni che puoi usare. Inoltre, offerte scikit-learn diversi metodi per eseguire una selezione euristica delle funzionalità, a condizione che il problema sia adatto alle loro tecniche. Ma trovare una buona euristica di scopo generale – nella forma più generale – è un problema difficile. In questa serie di articoli esploreremo alcune opzioni che potrebbero funzionare abbastanza bene anche quando N è piuttosto grande e la funzione obiettivo può essere letteralmente qualsiasi cosa, purché sia facile da calcolare.
Per tutte le tecniche di ottimizzazione presenti in questa serie di articoli, utilizzo il popolarissimo set di dati sui prezzi delle case su Kaggle (licenza MIT) — dopo alcune semplici trasformazioni di caratteristiche, finisce per avere 213 caratteristiche (N=213) e 1453 osservazioni. Il modello che sto utilizzando è la regressione lineare, statsmodels.api.OLS() e la funzione obiettivo che sto cercando di minimizzare è BIC: il criterio di informazione bayesiano, una misura della perdita di informazioni, quindi un BIC inferiore è migliore. È simile all’AIC – il criterio informativo di Akaike – tranne che il BIC tende a produrre modelli più parsimoniosi (preferisce modelli con meno caratteristiche). Ridurre al minimo l’AIC o il BIC tende a ridurre l’overfitting. Ma potrebbero essere utilizzate anche altre funzioni obiettivo, ad esempio R-quadrato (la varianza spiegata nell’obiettivo) o R-quadrato corretto: tieni solo presente che un R-quadrato più grande è migliore, quindi questo è un problema di massimizzazione.
In definitiva, la scelta della funzione obiettivo è qui irrilevante. Ciò che conta è che abbiamo una funzione obiettivo che cerchiamo costantemente di ottimizzare utilizzando varie tecniche.
Il codice completo utilizzato in questa serie di articoli è contenuto in un unico taccuino in il mio repository di selezione delle funzionalità – collegato anche alla fine. Fornirò qui frammenti di codice all’interno del testo, ma ti invitiamo a controllare il contesto completo nel taccuino.
Ma prima, esaminiamo una tecnica di selezione delle caratteristiche ben nota e collaudata, che utilizzeremo come confronto di base con le tecniche più complesse descritte in seguito.
SFS, la versione avanzata, è un algoritmo semplice. Inizia provando a scegliere una singola caratteristica e seleziona quella caratteristica che minimizza la funzione obiettivo. Una volta selezionata una funzione, rimane selezionata per sempre. Quindi prova ad aggiungervi un’altra funzionalità (per un totale di 2 funzionalità), in modo tale da minimizzare nuovamente l’obiettivo. Aumenta il numero di funzionalità selezionate di 1, cercando ogni volta di trovare la migliore nuova funzionalità da aggiungere alla raccolta esistente. La ricerca termina quando tutte le funzionalità vengono provate insieme. Qualunque sia la combinazione che minimizza l’obiettivo, vince.
SFS è un algoritmo avido – ogni scelta è ottimale a livello locale – e non torna mai indietro per correggere i propri errori. Ma è ragionevolmente veloce anche quando N è abbastanza grande. Il numero totale di combinazioni che tenta è N(N+1)/2 che è un polinomio quadratico (mentre una ricerca esaustiva richiede l’esecuzione di un numero esponenziale di prove).
Vediamo come può apparire il codice SFS in Python, utilizzando la libreria mlxtend:
import statsmodels.api as sm
from mlxtend.feature_selection import SequentialFeatureSelector as SFS
from sklearn.base import BaseEstimatorclass DummyEstimator(BaseEstimator):
# mlxtend wants to use an sklearn estimator, which is not needed here
# (statsmodels OLS is used instead)
# create a dummy estimator to pacify mlxtend
def fit(self, X, y=None, **kwargs):
return self
def neg_bic(m, X, y):
# objective function
lin_mod_res = sm.OLS(y, X, hasconst=True).fit()
return -lin_mod_res.bic
seq_selector = SFS(
DummyEstimator(),
k_features=(1, X.shape(1)),
forward=True,
floating=False,
scoring=neg_bic,
cv=None,
n_jobs=-1,
verbose=0,
# make sure the intercept is not dropped
fixed_features=('const'),
)
n_features = X.shape(1) - 1
objective_runs_sfs = round(n_features * (n_features + 1) / 2)
t_start_seq = time.time()
# mlxtend will mess with your dataframes if you don't .copy()
seq_res = seq_selector.fit(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()
Esegue rapidamente le combinazioni e questi sono i risultati riepilogativi:
best k: 36
best objective: 33708.98602877906
R2 @ best k: 0.9075677543649224
objective runs: 22791
total run time: 44.850 sec
Il numero migliore di funzioni è 36 su 213. Il BIC migliore è 33708.986 e richiede meno di 1 minuto per essere completato sul mio sistema. Invoca la funzione obiettivo 22,8k volte.
Questi sono i migliori valori BIC e R quadrato, in funzione del numero di caratteristiche selezionate:
Per ulteriori informazioni, come i nomi delle funzionalità effettivamente selezionate, controlla il notebook nel repository.
Ora proviamo qualcosa di più complesso.
Questo è un algoritmo di ottimizzazione numerica. Appartiene alla stessa classe degli algoritmi genetici (sono entrambi evolutivi), ma CMA-ES è abbastanza distinto da GA. È un algoritmo stocastico ed è privo di derivate: non richiede il calcolo di una derivata della funzione obiettivo (a differenza della discesa del gradiente, che si basa su derivate parziali). È efficiente dal punto di vista computazionale ed è utilizzato in una varietà di librerie di ottimizzazione numerica come Optuna. Tenterò qui solo un breve schizzo di CMA-ES; per spiegazioni più dettagliate si rimanda alla letteratura nella sezione link alla fine.
Considera la funzione Rastrigin bidimensionale:
La mappa termica seguente mostra i valori di questa funzione: colori più luminosi indicano un valore più alto. La funzione ha il minimo globale nell’origine (0, 0), ma è costellata di molti estremi locali. Vogliamo trovare il minimo globale tramite CMA-ES.
CMA-ES si basa sulla distribuzione normale multivariata. Genera punti di test nello spazio di ricerca da questa distribuzione. Dovrai indovinare il valore medio originale della distribuzione e la sua deviazione standard, ma dopodiché l’algoritmo modificherà iterativamente tutti questi parametri, spostando la distribuzione nello spazio di ricerca, cercando i migliori valori della funzione obiettivo. Ecco la distribuzione originale da cui vengono estratti i punti di test:
xi è l’insieme di punti che l’algoritmo genera ad ogni passaggio, nello spazio di ricerca. lambda è il numero di punti generati. Il valore medio della distribuzione verrà aggiornato ad ogni passo e, si spera, alla fine convergerà sulla vera soluzione. sigma è la deviazione standard della distribuzione: la diffusione dei punti di test. C è la matrice di covarianza: definisce la forma della distribuzione. A seconda dei valori di C la distribuzione può avere forma “tonda” oppure forma più allungata, ovale. Cambia in C consentire a CMA-ES di “intrufolarsi” in determinate aree dello spazio di ricerca o di evitare altre aree.
Nell’immagine sopra è stata generata una popolazione di 6 punti, che è la dimensione predefinita della popolazione scelta dall’ottimizzatore per questo problema. questo è il primo passo. Successivamente, l’algoritmo deve:
- calcolare la funzione obiettivo (Rastrigin) per ciascun punto
- aggiornare la media, la deviazione standard e la matrice di covarianza, creando di fatto una nuova distribuzione normale multivariata basata su ciò che ha imparato dalla funzione obiettivo
- generare una nuova serie di punti di test dalla nuova distribuzione
- ripetere finché non viene soddisfatto un criterio (convergere su un valore medio o superare il numero massimo di passaggi, ecc.)
Non mostrerò qui gli aggiornamenti per tutti i parametri di distribuzione, altrimenti questo articolo diventerà molto lungo: controlla i collegamenti alla fine per una spiegazione completa. Ma aggiornare solo la media della distribuzione è abbastanza semplice e funziona in questo modo: dopo aver calcolato la funzione obiettivo per ciascun punto del test, vengono assegnati i pesi ai punti, con i pesi maggiori assegnati ai punti con un valore obiettivo migliore e dalle loro posizioni viene calcolata una somma ponderata, che diventa la nuova media. In effetti, CMA-ES sposta la media della distribuzione verso i punti con un valore obiettivo migliore:
Se l’algoritmo converge alla soluzione vera, allora il valore medio della distribuzione convergerà a quella soluzione. La deviazione standard convergerà a 0. La matrice di covarianza cambierà la forma della distribuzione (rotonda o ovale), a seconda della geografia della funzione obiettivo, estendendosi in aree promettenti e allontanandosi da aree negative.
Ecco una GIF animata che mostra l’evoluzione nel tempo dei punti di prova mentre CMA-ES risolve il problema Rastrigin:
La funzione 2D Rastrigin era relativamente semplice, poiché ha solo 2 dimensioni. Per il nostro problema di selezione delle caratteristiche, abbiamo N=213 dimensioni. Inoltre lo spazio non è continuo. Ciascun punto di test è un vettore N-dimensionale con valori componenti da {0, 1} . In altre parole, ciascun punto di test si presenta così: (1, 0, 0, 1, 1, 1, 0, ...) – un vettore binario. Ma a parte questo, il problema è lo stesso: dobbiamo trovare i punti (i vettori) che minimizzano la funzione obiettivo: il parametro BIC di un modello OLS.
Ecco una versione semplice del codice CMA-ES per la selezione delle funzionalità, utilizzando la biblioteca del cmaes:
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).fit()
return lin_mod.bicX_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(
mean=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(range(max_gen)):
solutions = ()
for _ in range(optimizer.population_size):
x_for_eval, x_for_tell = optimizer.ask()
value = cma_objective(x_for_eval)
solutions.append((x_for_tell, value))
if value < best_objective_cmaes_small:
best_objective_cmaes_small = value
best_sol_raw_cmaes_small = x_for_eval
optimizer.tell(solutions)
best_features_cmaes_small = (
features_select(i)
for i, val in enumerate(best_sol_raw_cmaes_small.tolist())
if val == 1.0
)
print(f'best objective: {best_objective_cmaes_small}')
print(f'best features: {best_features_cmaes_small}')
L’ottimizzatore CMA-ES è definito con alcune ipotesi iniziali per il valore medio e per sigma (deviazione standard). Quindi esegue il ciclo attraverso molte generazioni, creando punti di test x_for_eval valutandoli con l’obiettivo, modificando la distribuzione (media, sigma, matrice di covarianza), ecc. Ciascuno x_for_eval il punto è un vettore binario (1, 1, 1, 0, 0, 1, ...) utilizzato per selezionare le caratteristiche dal set di dati.
Si prega di notare che CMAwM() viene utilizzato l’ottimizzatore (CMA con margine) invece del valore predefinito CMA(). L’ottimizzatore predefinito funziona bene per problemi regolari e continui, ma in questo caso lo spazio di ricerca è ad alta dimensione e sono consentiti solo due valori discreti (0 e 1). L’ottimizzatore predefinito rimane bloccato in questo spazio. CMAwM() allarga un po’ lo spazio di ricerca (mentre le soluzioni che restituisce sono ancora vettori binari), e questo sembra sufficiente per sbloccarlo.
Questo semplice codice funziona sicuramente, ma è tutt’altro che ottimale. Nel notebook associato ne ho una versione più complessa e ottimizzata, che è in grado di trovare soluzioni migliori, più velocemente. Ma il codice è grande, quindi non lo mostrerò qui — controlla il taccuino.
L’immagine seguente mostra la storia del codice CMA-ES complesso e ottimizzato alla ricerca della soluzione migliore. La mappa termica mostra la prevalenza/popolarità di ciascuna caratteristica in ciascuna generazione (più luminoso = più popolare). Puoi vedere come alcune funzionalità sono sempre molto apprezzate, mentre altre passano di moda rapidamente, e altre ancora vengono “scoperte” in seguito. La dimensione della popolazione scelta dall’ottimizzatore, dati i parametri di questo problema, è di 20 punti (individui), quindi la popolarità delle funzionalità viene calcolata in media su questi 20 punti.
Fonte: towardsdatascience.com
