AutoML

A esta altura, todo modelo que você treina vem com botões. Seleção de Modelos explicou o que os botões fazem (viés contra variância) e pôs várias famílias de modelos num único GridSearchCV; Validação Cruzada mostrou como comparar os resultados com honestidade e por que o vencedor de uma busca sempre parece melhor do que é. O que nenhuma das duas aulas perguntou é como gastar o orçamento da busca. Uma grade testa todas as combinações de alguns valores; com cinco hiperparâmetros e cinco valores cada, são 3.125 validações cruzadas, a maioria em regiões que já se mostravam ruins depois das primeiras doze.

A otimização de hiperparâmetros trata a própria busca como um problema de otimização; o AutoML leva a mesma ideia ao pipeline inteiro: pré-processamento, família de modelo, hiperparâmetros e o ensemble final. Esta aula constrói os métodos por dentro, cada um medido: busca aleatória contra uma grade (o argumento de Bergstra & Bengio, simulado), otimização bayesiana com um processo gaussiano e a melhoria esperada calculada à mão, TPE (tree-structured Parzen estimator) e Optuna, halving sucessivo e Hyperband. Depois, mede as duas maneiras como uma busca automática falha: ela otimiza um vazamento com o mesmo empenho que um sinal; quanto mais tempo busca, mais o próprio escore a lisonjeia.

Todo número abaixo vem dos blocos de código ou dos scripts das figuras que ficam ao lado desta página.

Como esta aula está organizada

A aula tem duas abordagens. Cada uma tem slides, texto e um lab para fazer em aula, no Jupyter ou no Colab; o quiz do fim cobre a fundamental. Comece pela fundamental: a avançada parte dela.

Os dois labs usam o mesmo notebook do lab, pronto para o Jupyter ou o Colab (Arquivo → Fazer upload de notebook): a parte fundamental vai do passo 1 ao 5; a avançada, do passo 6 ao 10.

Abordagem fundamental

Uma pergunta guia esta parte: como gastar um orçamento de busca e como saber o que a busca entregou de verdade? A grade e a busca aleatória não lembram nada; a otimização bayesiana aprende com as tentativas passadas onde olhar em seguida; o halving sucessivo corta cedo as configurações ruins. Cada método é calculado num problema pequeno o bastante para conferir à mão e depois medido. A parte termina com as duas maneiras como toda busca falha: ela otimiza um vazamento com o mesmo empenho que um sinal; quanto mais tempo busca, mais o próprio escore a lisonjeia.

Esta parte segue um caminho em cinco etapas, desenhado abaixo: 1. enquadrar a busca; 2. grade ou aleatória; 3. uma busca que aprende; 4. cortar tentativas ruins; 5. conferir o que ela entrega. Cada etapa leva à sua seção e, nos slides, uma faixa embaixo mostra em que etapa você está.

Slides da abordagem fundamental

Clique nos slides e use as setas ← → para navegar; F abre em tela cheia e Esc mostra todos os slides. Abrir os slides numa aba nova

1. Ajustar é um problema de otimização

Uma caixa-preta cara e ruidosa

Etapa 1. Chame de \(\lambda\) uma configuração (uma taxa de aprendizado, um tamanho de árvore, um tamanho de folha...) e de \(f(\lambda)\) o seu escore de validação cruzada. Ajustar é encontrar o \(\lambda\) com o maior \(f(\lambda)\). Escrito assim, parece com os problemas do gradiente descendente, mas três coisas o tornam diferente. Eis uma chamada a \(f\) no kc1, um conjunto de dados da NASA com 2.109 módulos de software descritos por 21 métricas de código (linhas de código, complexidade ciclomática, contagens de operadores...), em que o alvo diz se um defeito foi encontrado depois no módulo:

import numpy as np
from sklearn.datasets import fetch_openml
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import StratifiedKFold, cross_val_score, train_test_split

kc1 = fetch_openml(data_id=1067, as_frame=True)            # 2 109 módulos de um sistema de voo da NASA
X, y = kc1.data, (kc1.target == 'true').astype(int)         # 1 = um defeito foi encontrado depois
X_dev, X_test, y_dev, y_test = train_test_split(X, y, test_size=0.2, stratify=y, random_state=0)
print(X_dev.shape, f"defective: {y_dev.mean():.1%}")


def evaluate(params, seed=0):
    """Uma chamada à caixa-preta: o ROC AUC em CV 5-fold de um modelo de boosting com estes ajustes."""
    cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=seed)
    model = HistGradientBoostingClassifier(random_state=0, **params)
    return cross_val_score(model, X_dev, y_dev, cv=cv, scoring='roc_auc').mean()


for name, params in [('defaults', {}),
                     ('small, slow trees', {'learning_rate': 0.03, 'max_leaf_nodes': 8,
                                            'min_samples_leaf': 50})]:
    runs = [evaluate(params, seed) for seed in range(5)]    # mesmos ajustes, 5 embaralhamentos dos folds
    print(f"{name:17s} {np.round(runs, 3)}  mean {np.mean(runs):.3f}")

(1687, 21) defective: 15.5%
defaults          [0.76  0.776 0.767 0.769 0.769]  mean 0.768
small, slow trees [0.795 0.806 0.804 0.803 0.805]  mean 0.803

• Cara. Cada chamada envolve 5 ajustes de modelo. Aqui eles levam uma fração de segundo; num projeto real, uma chamada pode levar minutos ou horas, então o orçamento é o número de chamadas, não a aritmética em volta delas.
• Ruidosa. A mesma configuração, validada com 5 embaralhamentos diferentes dos folds, marca qualquer coisa entre 0,760 e 0,776. A diferença entre as duas configurações (0,768 contra 0,803) é real, mas uma diferença de 0,01 entre outras duas talvez não seja.
• Uma caixa-preta. Não existe gradiente do escore de CV em relação a min_samples_leaf: a função só pode ser avaliada. Vários botões são inteiros ou categorias; alguns só existem quando outro botão assume certo valor.

Um método de ajuste é uma regra para escolher o próximo \(\lambda\) a avaliar, dados os que já foram avaliados. Os métodos abaixo só diferem nessa regra.

Desenhando o espaço de busca

Antes de qualquer método rodar, alguém precisa dizer quais valores são permitidos. Duas decisões importam mais do que a escolha do método.

Escala. Uma taxa de aprendizado de 0,001 e outra de 0,01 são tão diferentes quanto 0,1 e 1: o que importa é a razão. Amostrar uniformemente entre 0,001 e 1 ignora isso:

from scipy.stats import loguniform, uniform

u = uniform(0.001, 1 - 0.001).rvs(100_000, random_state=0)   # uniforme em [0.001, 1]
lu = loguniform(0.001, 1).rvs(100_000, random_state=0)       # uniforme na escala log
for lo, hi in [(0.001, 0.01), (0.01, 0.1), (0.1, 1)]:
    print(f"learning rate in [{lo}, {hi}):  uniform {np.mean((u >= lo) & (u < hi)):5.1%}"
          f"   log-uniform {np.mean((lu >= lo) & (lu < hi)):5.1%}")

learning rate in [0.001, 0.01):  uniform  0.9%   log-uniform 33.5%
learning rate in [0.01, 0.1):  uniform  9.1%   log-uniform 33.2%
learning rate in [0.1, 1):  uniform 90.0%   log-uniform 33.3%

À mão: o intervalo \([0{,}001;\ 0{,}01)\) é \(0{,}009 / 0{,}999 = 0{,}9\%\) da faixa uniforme, enquanto na escala log cada década é um terço das \(\log_{10}(1/0{,}001) = 3\) décadas. Uma busca uniforme gasta 90% do orçamento em taxas de aprendizado acima de 0,1 e quase nada na década em que, no kc1, vivem os melhores modelos de boosting (veja a figura do halving mais abaixo: as melhores configurações têm taxas de aprendizado entre 0,01 e 0,03). Regra: tudo o que age de forma multiplicativa (taxas de aprendizado, forças de regularização C e alpha, o número de folhas ou de vizinhos) ganha escala log: loguniform no scikit-learn, log=True no Optuna.

Parâmetros condicionais. Alguns hiperparâmetros só existem em algumas configurações: C pertence à regressão logística, max_leaf_nodes ao boosting, gamma a um kernel RBF. Uma grade única sobre todos eles combinaria cada C com cada tamanho de árvore, ajustando a mesma regressão logística repetidamente. Seleção de Modelos resolveu isso com uma lista de grades, uma por família. O espaço de busca é, na verdade, uma árvore: primeiro escolha a família, depois os botões dessa família. O código do Optuna da abordagem avançada escreve essa árvore como um if.

2. Grade ou aleatória? Bergstra & Bengio, medido

Etapa 2. Seleção de Modelos terminou com uma afirmação de Bergstra & Bengio (2012): quando só alguns hiperparâmetros importam de verdade, combinações aleatórias cobrem esses poucos muito melhor do que uma grade. Eis o argumento, medido. Tome um escore que depende de um hiperparâmetro, \(x_1\), por meio de uma lombada de largura 0,05 em torno de um melhor valor desconhecido \(c\); os outros \(d - 1\) hiperparâmetros não fazem nada. Os dois métodos recebem o mesmo orçamento \(B\). Uma grade com \(k\) valores por eixo custa \(B = k^d\) avaliações, mas testa o parâmetro importante em apenas \(k\) valores distintos. A busca aleatória o testa em \(B\) valores distintos, qualquer que seja \(d\).

À mão: com \(k\) valores de grade em \(\tfrac{1}{2k}, \tfrac{3}{2k}, \dots\), o ótimo fica no máximo a \(\tfrac{1}{2k}\) do valor mais próximo e, em média, a \(\tfrac{1}{4k}\): 0,083 para \(k = 3\). Com \(B\) valores aleatórios, o mais próximo fica em média a cerca de \(\tfrac{1}{2(B + 1)}\): 0,05 para \(B = 9\), 0,006 para \(B = 81\). Na média de 4.000 problemas com \(c\) sorteado:

rng = np.random.default_rng(0)
c = rng.uniform(0, 1, 4000)         # 4 000 problemas: o melhor valor de x1 pode estar em qualquer lugar
W = 0.05                            # o escore é uma lombada de largura 0.05 em torno dele


def distance_and_score(x1):
    """Distância média do ótimo ao x1 testado mais próximo e melhor escore médio encontrado."""
    dist = np.abs(x1 - c[:, None]).min(axis=1)
    return dist.mean(), np.exp(-dist ** 2 / (2 * W ** 2)).mean()


print(" d  budget  x1 values in grid   distance: grid  random   best score: grid  random")
for d, k in [(2, 3), (2, 9), (3, 3), (4, 3), (4, 4)]:
    B = k ** d                                          # uma grade k × k × ...
    grid_x1 = (np.arange(k) + 0.5) / k                  # os k valores que a grade dá a x1
    rand_x1 = rng.uniform(0, 1, (len(c), B))            # B valores aleatórios de x1, novos por problema
    (dg, sg), (dr, sr) = distance_and_score(grid_x1), distance_and_score(rand_x1)
    print(f"{d:2d}  {B:6d}  {k:17d}   {dg:14.3f}  {dr:6.3f}   {sg:16.3f}  {sr:6.3f}")

 d  budget  x1 values in grid   distance: grid  random   best score: grid  random
 2       9                  3            0.083   0.054              0.378   0.629
 2      81                  9            0.028   0.006              0.824   0.985
 3      27                  3            0.083   0.018              0.378   0.900
 4      81                  3            0.083   0.006              0.378   0.986
 4     256                  4            0.062   0.002              0.501   0.998

(script da figura)

Com 9 avaliações em duas dimensões, a grade cai em média a 0,083 do ótimo e marca 0,378; a busca aleatória cai a 0,054 e marca 0,629. A diferença explode com o número de dimensões inúteis. Com 81 avaliações sobre quatro hiperparâmetros, a grade ainda testa \(x_1\) em apenas 3 valores, então marca exatamente o que a grade 3 × 3 marcava (0,378), enquanto a busca aleatória está em 0,986. A grade paga cada eixo irrelevante com uma multiplicação do seu custo; a busca aleatória não paga nada, porque um eixo irrelevante não muda onde seus pontos caem no eixo relevante.

Buscas reais são assim. Bergstra & Bengio ajustaram um modelo do escore de validação aos próprios experimentos com redes neurais e descobriram que só alguns hiperparâmetros importavam de verdade em cada conjunto de dados e que quais importavam mudava de um conjunto para outro; por isso não dá para desenhar a grade em torno deles com antecedência. Uma grade continua boa para um ou dois hiperparâmetros que você entende (uma curva de C, uma curva de k). Além disso, a busca aleatória é o baseline honesto que todo método mais esperto precisa vencer.

Veja um problema de cada vez. Comece com orçamento 9, conte os valores distintos de \(x_1\) que cada método testa; depois aumente o orçamento e aperte novo problema algumas vezes. A linha de baixo, à direita, refaz a média sobre 2.000 problemas:

Num problema isolado a grade às vezes ganha, quando uma de suas poucas colunas cai por acaso no pico; na média ela nunca ganha, porque uma grade \(k \times k\) paga \(k^2\) avaliações por \(k\) valores do único hiperparâmetro que importa. O segundo modo do simulador, halving sucessivo, antecipa o torneio da seção 4: os escores são sintéticos, mas reproduzem a armadilha descrita lá, em que um orçamento contado em árvores elimina primeiro quem aprende devagar.

3. Otimização bayesiana: uma busca que aprende

Etapa 3. A busca aleatória ignora tudo o que já viu: a 27.ª configuração é sorteada exatamente como a primeira. Quando cada avaliação custa uma hora, isso é desperdício. A otimização bayesiana (Jones, Schonlau & Welch, 1998; Snoek, Larochelle & Adams, 2012, para aprendizado de máquina) usa as avaliações feitas até agora para decidir onde olhar em seguida, num laço:

1. ajuste um modelo substituto (surrogate) barato aos pares (configuração, escore) vistos até agora. O substituto prevê o escore em qualquer ponto, com uma incerteza;
2. escolha a próxima configuração maximizando uma função de aquisição, que pondera exploração local (exploitation: onde o substituto prevê escore alto) contra exploração (exploration: onde ele não tem certeza);
3. avalie essa configuração, acrescente ao histórico, repita.

Para vê-la funcionar, precisamos de uma caixa-preta cuja verdade conhecemos. A função abaixo faz o papel de um escore de validação em função de um hiperparâmetro \(x\) (em escala log): um pico largo de 0,78 em 0,25 e outro, mais alto e mais estreito, de 0,82 em 0,72.

O modelo substituto

O substituto usual é um processo gaussiano (GP). Você não precisa da matemática dele para usá-lo: um GP prevê, em cada \(x\), uma média \(\mu(x)\) e um desvio-padrão \(\sigma(x)\). A média passa pelos pontos avaliados; o desvio-padrão é zero neles e cresce com a distância a eles. A rapidez com que cresce é dada pela escala de comprimento (length scale) do kernel: com 0,15, o substituto supõe que pontos a menos de uns 0,15 um do outro têm escores parecidos.

from scipy.stats import norm
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF


def score(x):
    """Um substituto para um escore de validação caro, em função de um hiperparâmetro."""
    return 0.70 + 0.08 * np.exp(-((x - 0.25) / 0.10) ** 2) + 0.12 * np.exp(-((x - 0.72) / 0.06) ** 2)


X_seen = np.array([0.15, 0.35, 0.95])                       # três avaliações até agora
y_seen = score(X_seen)
gp = GaussianProcessRegressor(kernel=RBF(length_scale=0.15), normalize_y=True, optimizer=None)
gp.fit(X_seen.reshape(-1, 1), y_seen)
print('evaluated:', X_seen, '->', y_seen.round(5))
for x in [0.25, 0.55, 0.72]:
    mu, sd = gp.predict([[x]], return_std=True)
    print(f"x = {x}: surrogate mean {mu[0]:.5f}, sd {sd[0]:.5f}")

evaluated: [0.15 0.35 0.95] -> [0.72943 0.72943 0.7    ]
x = 0.25: surrogate mean 0.73076, sd 0.00419
x = 0.55: surrogate mean 0.72212, sd 0.01246
x = 0.72: surrogate mean 0.71390, sd 0.01318

Três pontos parecem muito diferentes para o substituto. Em 0,25, entre dois pontos avaliados, ele está quase certo: média 0,73076, logo acima do melhor escore até agora, \(f^+ = 0{,}72943\), com um desvio-padrão minúsculo. Em 0,55 e 0,72, longe de qualquer avaliação, as médias estão abaixo de \(f^+\), mas a incerteza é três vezes maior. (O substituto não sabe nada do pico em 0,72: nada foi avaliado perto dele.)

A função de aquisição: melhoria esperada, à mão

Qual dos três merece a próxima avaliação? A melhoria esperada (expected improvement, EI) pergunta: se o escore em \(x\) é uma variável normal com média \(\mu\) e desvio-padrão \(\sigma\), quanto ele vai superar \(f^+\) em média, contando como zero quando não supera? A resposta tem forma fechada, com \(z = (\mu - f^+)/\sigma\), \(\Phi\) a distribuição acumulada normal e \(\varphi\) a sua densidade:

O primeiro termo premia uma média alta (exploração local), o segundo uma incerteza grande (exploração). Com os números impressos acima (\(\Phi\) e \(\varphi\) de uma tabela da normal):

Em 0,25 quem trabalha é a média; em 0,55 é a incerteza, que quase vence. A biblioteca concorda:

def expected_improvement(x, gp, best):
    mu, sd = gp.predict(np.reshape(x, (-1, 1)), return_std=True)
    z = (mu - best) / sd
    return (mu - best) * norm.cdf(z) + sd * norm.pdf(z)


xs = np.linspace(0, 1, 201)                                 # os candidatos que podemos pontuar
ei = expected_improvement(xs, gp, y_seen.max())
for x in [0.25, 0.55, 0.72]:
    print(f"EI({x}) = {expected_improvement(x, gp, y_seen.max())[0]:.4f}")
print(f"largest EI at x = {xs[ei.argmax()]:.3f}")

EI(0.25) = 0.0024
EI(0.55) = 0.0021
EI(0.72) = 0.0008
largest EI at x = 0.250

Calcular a EI em 201 candidatos não custa nada perto de uma validação cruzada: essa é toda a troca. A função cara é avaliada uma vez por iteração; o substituto barato é consultado milhares de vezes.

Quatro iterações

def bayes_opt(X_seen, n_iter, length_scale=0.15):
    """Ajusta o substituto, avalia o argmax da EI, repete."""
    y_seen = score(X_seen)
    for _ in range(n_iter):
        gp = GaussianProcessRegressor(kernel=RBF(length_scale), normalize_y=True, optimizer=None)
        gp.fit(X_seen.reshape(-1, 1), y_seen)
        x_next = xs[np.argmax(expected_improvement(xs, gp, y_seen.max()))]
        X_seen, y_seen = np.append(X_seen, x_next), np.append(y_seen, score(x_next))
    return X_seen, y_seen


X_bo, y_bo = bayes_opt(np.array([0.15, 0.35, 0.95]), n_iter=4)
print('x tried:', X_bo[3:], '-> score', y_bo[3:].round(3))

x tried: [0.25  0.72  0.685 0.775] -> score [0.78  0.82  0.785 0.752]

(script da figura)

Cada coluna é uma iteração: o substituto (faixa azul: média ± 2 desvios-padrão) com o escore verdadeiro em cinza tracejado e, embaixo, a melhoria esperada com o máximo marcado. A iteração 1 explora localmente: avalia 0,25 e encontra 0,780. Agora o pico da esquerda é conhecido, então a EI ali desaba e a maior EI passa para o vão largo e inexplorado entre 0,4 e 0,9. A iteração 2 explora ali, avalia 0,72 e encontra o ótimo, 0,820. As iterações 3 e 4 conferem os vizinhos, 0,685 e 0,775. Com cinco avaliações ao todo, a busca encontrou a melhor configuração; a busca aleatória precisaria, em média, de uns 28 sorteios para cair a menos de \(\pm 0{,}018\) de 0,72, a faixa em que o escore passa de 0,81.

O simulador refaz as quatro iterações, um passo de cada vez: em cada uma, a média e a incerteza do substituto nas candidatas, a melhoria esperada escrita com esses números, o ponto escolhido e o escore que ele devolve. Antes de apertar próximo, preveja para onde vai a iteração 2 depois que 0,25 for avaliado.

Aquela rodada foi boa. A abordagem avançada a repete a partir de 100 partidas aleatórias e mostra como as hipóteses do próprio substituto decidem o que ele consegue achar.

4. Halving sucessivo: cortar cedo as tentativas ruins

Etapa 4. A otimização bayesiana gasta bem o orçamento escolhendo onde olhar. O halving sucessivo economiza de outro jeito: olha muitas configurações, mas gasta pouco com cada uma até ter evidência de que ela vale mais.

Um torneio, à mão

A maioria das configurações ruins se mostra ruim cedo: depois de poucas rodadas de boosting ou treinada numa amostra pequena. O halving sucessivo (successive halving; Jamieson & Talwalkar, 2016) organiza um torneio. Comece \(n\) configurações com um orçamento pequeno \(r\); fique com a melhor fração \(1/\eta\) delas; dê aos sobreviventes \(\eta\) vezes mais orçamento; repita até sobrar uma. Com \(n = 27\), \(\eta = 3\) e unidades de orçamento de 10 árvores:

Toda rodada custa o mesmo, 270 ajustes de árvore, então o torneio inteiro custa \(4 \times 270 = 1\,080\). Treinar as 27 configurações com 270 árvores custaria \(27 \times 270 = 7\,290\): o torneio olha 27 configurações pelo preço de 4 completas. O scikit-learn o implementa como HalvingRandomSearchCV (ainda experimental, daí o import extra). A unidade de orçamento pode ser o número de linhas (o padrão) ou qualquer hiperparâmetro inteiro, como o número de árvores:

import pandas as pd
from scipy.stats import randint
from sklearn.experimental import enable_halving_search_cv  # noqa: F401  (ainda experimental)
from sklearn.model_selection import HalvingRandomSearchCV, ParameterSampler

space = {'learning_rate': loguniform(0.01, 0.5), 'max_leaf_nodes': randint(4, 64),
         'min_samples_leaf': randint(5, 200), 'l2_regularization': loguniform(1e-3, 10)}
base = HistGradientBoostingClassifier(max_iter=270, random_state=0)
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)       # um só cv para todas as buscas abaixo
halving = {}
for budget, kw in [('trees', dict(resource='max_iter', min_resources=10, max_resources=270)),
                   ('rows', dict(min_resources=60))]:
    halving[budget] = HalvingRandomSearchCV(base, space, n_candidates=27, factor=3, cv=cv,
                                            scoring='roc_auc', random_state=0, **kw).fit(X_dev, y_dev)
    print(f"budget = {budget:5s}: candidates {halving[budget].n_candidates_}, "
          f"resources {halving[budget].n_resources_}")

budget = trees: candidates [27, 9, 3, 1], resources [10, 30, 90, 270]
budget = rows : candidates [27, 9, 3, 1], resources [60, 180, 540, 1620]

A biblioteca rodou exatamente a tabela: 27, 9, 3 e 1 candidatos, com 10, 30, 90 e 270 árvores. A segunda rodada é o mesmo torneio com linhas como orçamento: 60, 180, 540 e 1.620 das 1.687 linhas de desenvolvimento.

O que o torneio supõe

O halving só é tão bom quanto a hipótese central dele: a ordem depois de um orçamento pequeno prevê a ordem depois do orçamento completo. Para conferir, pontue os 27 candidatos com o orçamento completo (270 árvores, todas as linhas) e veja onde cai o vencedor de cada torneio:

from sklearn.base import clone

full = {}                                               # todo candidato com o orçamento completo
for p in ParameterSampler(space, n_iter=27, random_state=0):   # os 27 que o halving sorteou
    full[p['learning_rate']] = cross_val_score(clone(base).set_params(**p), X_dev, y_dev, cv=cv,
                                               scoring='roc_auc').mean()
full = pd.Series(full).sort_values(ascending=False)
print(f"best at the full budget: {full.iloc[0]:.3f} (learning rate {full.index[0]:.3f})")
for budget, h in halving.items():
    lr = h.best_params_['learning_rate']
    print(f"winner with budget = {budget:5s}: learning rate {lr:.3f}, full-budget score {full[lr]:.3f}, "
          f"rank {full.index.get_loc(lr) + 1} of 27")

best at the full budget: 0.791 (learning rate 0.016)
winner with budget = trees: learning rate 0.062, full-budget score 0.781, rank 13 of 27
winner with budget = rows : learning rate 0.013, full-budget score 0.787, rank 5 of 27

(script da figura)

Com linhas como orçamento, o vencedor fica em 5.º de 27 com o orçamento completo, 0,787 contra o melhor 0,791: uma boa escolha por um sétimo do custo. Com árvores como orçamento, o vencedor fica em 13.º (0,781); a figura mostra por quê. As melhores configurações com 270 árvores têm todas taxas de aprendizado pequenas (0,01–0,03) e todas foram eliminadas na primeira rodada. Você viu o motivo em Gradient Boosting: com uma taxa de aprendizado pequena, 10 árvores mal começaram, então quem aprende devagar sempre parece pior no começo do que quem aprende rápido, mesmo quando termina melhor. Quando a unidade de orçamento interage com um hiperparâmetro, a ordem inicial é enviesada e o halving elimina as configurações que deveria manter. Duas saídas: usar uma unidade de orçamento que não interaja com os botões buscados (linhas se saíram bem melhor aqui, embora, como mostra o simulador abaixo, não perfeitamente) ou usar parada antecipada e deixar a taxa de aprendizado e o número de árvores serem ajustados juntos.

O simulador joga os dois torneios rodada por rodada, nas mesmas 27 configurações: o escore de cada sobrevivente com o orçamento da rodada, o terço que segue, o custo de cada rodada e, no fim, onde cada configuração fica de verdade com o orçamento completo. Antes de avançar, preveja quais taxas de aprendizado sobrevivem à primeira rodada quando o orçamento é contado em árvores.

O torneio em linhas tem uma fraqueza própria, visível nas duas primeiras rodadas. Com 60 linhas (48 por fold de treino), 26 das 27 configurações pedem folhas de mais de 24 linhas: nem uma divisão é possível; as 26 marcam exatamente 0,500. Com 180 linhas, o mesmo acontece com 7 dos 9 sobreviventes. Esses cortes são decididos pelo jeito de desempatar, não pelos dados; o vencedor final (taxa de aprendizado 0,013, folhas de pelo menos 179 linhas) passou pelos dois num empate: uma boa escolha e em parte sortuda. O tamanho da folha interage com o número de linhas assim como a taxa de aprendizado interage com o número de árvores. O Hyperband, na abordagem avançada, divide essa aposta rodando vários torneios que começam com orçamentos diferentes.

5. O que o AutoML não consegue fazer

Etapa 5. Uma busca automática otimiza exatamente o número que você lhe dá, exatamente nos dados que você lhe dá. Tudo o que estiver errado com esse número ou com esses dados, ela otimiza também — e com mais força do que uma pessoa faria.

Ele otimiza o vazamento que você deu

Suponha que o kc1 viesse com mais uma coluna: bug_reports, o número de relatos de bug abertos contra cada módulo. É uma medida real e fortemente ligada ao alvo. Mas ela é preenchida depois do lançamento, enquanto o modelo existe para apontar módulos arriscados antes do lançamento, quando a contagem ainda é zero para todo mundo. É o atributo que conhece o futuro. Dê à mesma busca aleatória as duas versões da tabela:

from sklearn.metrics import roc_auc_score
from sklearn.model_selection import RandomizedSearchCV

rng = np.random.default_rng(0)


def add_bug_reports(X, y):
    """Uma coluna preenchida DEPOIS do lançamento: relatos de bug abertos contra o módulo."""
    X = X.copy()
    X['bug_reports'] = rng.poisson(np.where(y == 1, 2.0, 0.3))
    return X


X_dev_leak, X_test_leak = add_bug_reports(X_dev, y_dev), add_bug_reports(X_test, y_test)
X_test_prod = X_test_leak.assign(bug_reports=0)     # na hora da previsão, nada foi aberto ainda

searches = {}
for name, X_search in [('without the column', X_dev), ('with bug_reports', X_dev_leak)]:
    searches[name] = RandomizedSearchCV(HistGradientBoostingClassifier(random_state=0), space,
                                        n_iter=20, cv=cv, scoring='roc_auc',
                                        random_state=0).fit(X_search, y_dev)
    print(f"{name:19s} best CV ROC AUC {searches[name].best_score_:.3f}")


def test_auc(search, X_):
    return roc_auc_score(y_test, search.predict_proba(X_)[:, 1])


print(f"test, honest winner:                     {test_auc(searches['without the column'], X_test):.3f}")
print(f"test, leaky winner, column filled in:    {test_auc(searches['with bug_reports'], X_test_leak):.3f}")
print(f"test, leaky winner, as it will be used:  {test_auc(searches['with bug_reports'], X_test_prod):.3f}")

without the column  best CV ROC AUC 0.792
with bug_reports    best CV ROC AUC 0.916
test, honest winner:                     0.804
test, leaky winner, column filled in:    0.962
test, leaky winner, as it will be used:  0.782

A busca faz o trabalho dela perfeitamente. Com a coluna, a melhor configuração chega a um ROC AUC de CV de 0,916 (contra 0,792) e a 0,962 num conjunto de teste montado do mesmo jeito. Toda verificação automática passa, porque toda verificação rodou na tabela histórica, em que a coluna está preenchida. Nos módulos como eles serão de fato pontuados, com bug_reports = 0, o mesmo modelo cai para 0,782, abaixo dos 0,804 do modelo honesto. Ele aprendeu a se apoiar numa coluna que não vai existir e ficou pior do que um modelo que nunca a viu. Nenhum método de busca, por mais esperto, pega isso: o vazamento está nos dados. A importância por permutação do modelo final mostraria em que coluna ele se apoia, mas só alguém que sabe quando cada coluna é registrada pode dizer que se apoiar nela é ilegítimo.

Uma busca mais longa sobreajusta os dados de validação

A maldição do vencedor diz que o melhor de muitos escores ruidosos é otimista. Uma rodada de AutoML é uma máquina de produzir muitos escores ruidosos, então a maldição cresce com o orçamento. A figura roda muitas vezes uma busca aleatória de até 150 tentativas no kc1 e acompanha a configuração que ela escolheria: o escore de CV dela (o que a busca relata) e o escore dela num conjunto de teste grande e intocado (o que ela entrega) (script da figura):

Com 1.687 linhas de desenvolvimento, vale a pena buscar por mais tempo, como a busca informa com honestidade: de 10 para 150 tentativas, o escore de CV da configuração escolhida vai de 0,803 para 0,808 e o escore de teste, de 0,804 para 0,809. (A curva de teste começa acima da de CV porque o modelo final é reajustado em todas as linhas de desenvolvimento, 25% a mais do que cada fold de CV usou para treinar.) Com 300 linhas de desenvolvimento, a história é outra. As primeiras 10 tentativas ajudam muito (a primeira configuração aleatória costuma ser ruim), mas de 10 para 150 tentativas o escore relatado sobe 0,016, de 0,772 para 0,788, enquanto o entregue se move 0,004, de 0,765 para 0,769. Três quartos do ganho relatado são a busca se ajustando às peculiaridades de 300 linhas. Quanto menos dados, mais ruidoso é cada escore de CV e mais uma busca longa encontra configurações sortudas em vez de boas. Duas consequências:

• Orçamento não é de graça. A partir de certo ponto, mais tentativas compram otimismo, não desempenho. Pare quando o melhor escore deixar de se mover mais do que o próprio desvio-padrão entre folds.
• Avalie a busca inteira, não o vencedor dela. O escore que uma ferramenta de AutoML relata para o melhor modelo encontrado carrega a maldição. Meça-o num conjunto de teste trancado, uma vez, ou envolva a busca inteira numa validação cruzada externa (CV aninhada), como faz o lab avançado.

6. Resumo: ajuste em cinco linhas

1. Ajustar é uma caixa-preta cara e ruidosa. O orçamento é o número de chamadas; uma diferença menor que o ruído entre folds pode não ser real.
2. Desenhe o espaço antes do método. Escala log para tudo o que é multiplicativo; um botão pertence à sua família, então o espaço é uma árvore.
3. A busca aleatória vence a grade quando poucos botões importam: a grade paga cada eixo irrelevante, a busca aleatória não. A otimização bayesiana aprende: um substituto prevê o escore com uma incerteza e a melhoria esperada pondera um escore previsto alto contra uma incerteza grande.
4. O halving sucessivo olha muitas configurações pelo preço de poucas, se a ordem com pouco orçamento prevê a ordem com o orçamento completo; uma unidade de orçamento que interage com um botão (árvores e a taxa de aprendizado) quebra isso.
5. Uma busca otimiza exatamente o que você lhe dá: um vazamento com o mesmo empenho que um sinal e, quanto mais roda, o ruído dos próprios folds. Avalie-a em dados que ela nunca tocou.

Lab fundamental — grade ou aleatória, nas vogais

Os dados são o phoneme, de um projeto de reconhecimento de fala: 5.404 sons de vogais, cada um descrito pelas amplitudes dos cinco primeiros harmônicos. A tarefa é distinguir vogais orais de nasais. O modelo é um HistGradientBoostingClassifier com parada antecipada; a métrica é o ROC AUC. Esta parte carrega os dados, tranca um conjunto de teste, define as baselines que toda busca precisa bater, escreve o espaço de busca e roda as duas buscas que não lembram nada, uma grade e a busca aleatória, com o mesmo orçamento: 27 avaliações de uma validação cruzada 3-fold cada. O lab avançado acrescenta o TPE e o halving sucessivo, testa o vencedor uma vez e avalia a própria busca.

Passo 1 — carregar as vogais

Preveja. Vogais orais e nasais: você espera as duas classes em números parecidos? E que ROC AUC tem um modelo que ordena as vogais ao acaso?

# No Colab, instale o Optuna primeiro:  !pip install optuna
import numpy as np
import pandas as pd
from sklearn.datasets import fetch_openml

phoneme = fetch_openml(data_id=1489, as_frame=True)    # 5 404 vogais, 5 amplitudes de harmônicos cada
X = phoneme.data
y = (phoneme.target == '2').astype(int)                 # 1 = vogal oral, 0 = vogal nasal
print(X.shape)
print(y.value_counts())
print(X.describe().loc[['mean', 'std', 'min', 'max']].round(2))

(5404, 5)
Class
0    3818
1    1586
Name: count, dtype: int64
        V1    V2    V3    V4    V5
mean  0.00  0.00  0.00 -0.00 -0.00
std   1.00  1.00  1.00  1.00  1.00
min  -2.93 -3.04 -2.80 -2.49 -2.37
max   3.83  3.66  2.63  3.05  4.59

5.404 vogais: 1.586 orais (29%) e 3.818 nasais. As cinco colunas são amplitudes de harmônicos que os autores do conjunto já padronizaram (média 0, desvio-padrão 1); de qualquer forma, um modelo de boosting não precisaria de escala, então não há passo de pré-processamento.

Passo 2 — trancar o teste, fixar os folds, definir baselines

Preveja. Uma regressão logística traça uma única fronteira reta no espaço das cinco amplitudes; um modelo de boosting soma muitas caixinhas. Qual você espera mais alto? Por muito ou por pouco? Que diferença entre duas configurações você ainda chamaria de ruído?

from sklearn.dummy import DummyClassifier
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold, cross_val_score, train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X_dev, X_test, y_dev, y_test = train_test_split(X, y, test_size=0.2, stratify=y, random_state=42)
cv = StratifiedKFold(n_splits=3, shuffle=True, random_state=42)     # um único cv para todas as buscas


def model(**params):
    """O modelo que toda busca ajusta; a parada antecipada escolhe o número de árvores."""
    return HistGradientBoostingClassifier(max_iter=500, early_stopping=True, n_iter_no_change=10,
                                          random_state=42, **params)


baselines = {'dummy': DummyClassifier(),
             'logistic regression': make_pipeline(StandardScaler(), LogisticRegression()),
             'boosting, defaults': model()}
baseline = {}
for name, est in baselines.items():
    s = cross_val_score(est, X_dev, y_dev, cv=cv, scoring='roc_auc')
    baseline[name] = s.mean()
    print(f"{name:20s} ROC AUC {s.mean():.4f} ± {s.std():.4f}")

dummy                ROC AUC 0.5000 ± 0.0000
logistic regression  ROC AUC 0.8195 ± 0.0154
boosting, defaults   ROC AUC 0.9468 ± 0.0042

O dummy fica em 0,5, como sempre no ROC AUC. A regressão logística chega a 0,8195: nessas cinco amplitudes, as duas classes estão longe de ser linearmente separáveis. O modelo de boosting com os valores padrão (e parada antecipada) chega a 0,9468 ± 0,0042: esse é o número que toda busca precisa vencer. Guarde o desvio-padrão: diferenças de uns 0,004 entre configurações são do tamanho do ruído de fold para fold. model() é o único estimador que as buscas vão ajustar e cv é compartilhado por todas elas; assim, toda configuração de todo método é pontuada nos mesmos três folds.

Passo 3 — o espaço de busca

Preveja. A taxa de aprendizado é sorteada em escala log entre 0,01 e 0,5. A maioria dos sorteios será grande, acima de 0,1, como daria um sorteio uniforme, ou vai se espalhar pelas décadas?

from scipy.stats import loguniform, rv_discrete
from sklearn.model_selection import ParameterSampler


def log_int(low, high):
    """Inteiros de low a high, uniformes na escala log: P(k) proporcional a 1/k."""
    k = np.arange(low, high + 1)
    return rv_discrete(values=(k, (1 / k) / (1 / k).sum()))


space = {'learning_rate': loguniform(0.01, 0.5),
         'max_leaf_nodes': log_int(4, 128),
         'min_samples_leaf': log_int(2, 100),
         'l2_regularization': loguniform(1e-3, 10),
         'max_features': [0.5, 0.75, 1.0]}
BUDGET = 27                                             # configurações por método, 3 ajustes cada


def show(params):
    return ', '.join(f"{k}={v:.3g}" for k, v in sorted(params.items()))


for p in ParameterSampler(space, n_iter=4, random_state=0):
    print(show(p))

l2_regularization=0.157, learning_rate=0.164, max_features=0.75, max_leaf_nodes=33, min_samples_leaf=8
l2_regularization=0.0155, learning_rate=0.0125, max_features=0.5, max_leaf_nodes=14, min_samples_leaf=42
l2_regularization=0.13, learning_rate=0.0923, max_features=0.75, max_leaf_nodes=71, min_samples_leaf=6
l2_regularization=0.391, learning_rate=0.0422, max_features=0.75, max_leaf_nodes=58, min_samples_leaf=58

Quatro sorteios do espaço. rv_discrete monta uma distribuição de inteiros a partir de uma lista de valores e das suas probabilidades; com probabilidades proporcionais a \(1/k\), a faixa de 4–8 folhas é tão provável quanto a de 64–128, a escala log de Desenhando o espaço de busca para inteiros. Qualquer objeto com um método .rvs pode entrar num espaço de busca do scikit-learn. As faixas são largas de propósito: um fator de 50 na taxa de aprendizado, de 4 a 128 folhas, de 2 a 100 linhas por folha. O número de árvores não está no espaço: a parada antecipada o escolhe. Cada método recebe BUDGET = 27 configurações, isto é, 81 ajustes do modelo.

Passo 4 — busca em grade

Preveja. Uma grade 3 × 3 × 3 nos três botões que todo mundo ajusta. Você espera o vencedor no meio da grade ou numa borda? O que uma borda lhe diria?

from sklearn.model_selection import GridSearchCV

grid = {'learning_rate': [0.03, 0.1, 0.3],
        'max_leaf_nodes': [8, 32, 128],
        'min_samples_leaf': [5, 20, 80]}                # 3 × 3 × 3 = 27 configurações
grid_search = GridSearchCV(model(), grid, cv=cv, scoring='roc_auc').fit(X_dev, y_dev)
print(f"grid:   best ROC AUC {grid_search.best_score_:.4f} · {grid_search.best_params_}")

grid:   best ROC AUC 0.9498 · {'learning_rate': 0.3, 'max_leaf_nodes': 128, 'min_samples_leaf': 5}

Uma grade razoável, do tipo que qualquer um escreveria à mão: três valores dos três botões que todo mundo ajusta. O vencedor dela chega a 0,9498 e fica num canto da grade: a maior taxa de aprendizado, o maior número de folhas e o menor tamanho de folha. Um vencedor na borda de uma grade diz que a grade era estreita demais naquela direção; para descobrir o que há além, seria preciso estendê-la e pagar outra rodada de avaliações.

Simulação: uma grade montada em volta do botão que importa. A grade do passo 4 gastou um terço do orçamento em árvores de 8 folhas e achou o vencedor no canto das árvores grandes e folhas pequenas. Gaste as mesmas 27 avaliações numa grade montada em volta dessa descoberta: 9 valores de max_leaf_nodes, de 16 a 256, as mesmas 3 taxas de aprendizado e min_samples_leaf fixo em 5. Ela bate 0,9498? Qual número de folhas vence? Você conseguiria ter escrito essa grade antes de rodar o passo 4?

leaves_grid = {'learning_rate': [0.03, 0.1, 0.3],
               'max_leaf_nodes': [16, 24, 32, 48, 64, 96, 128, 192, 256],
               'min_samples_leaf': [5]}                    # 3 × 9 × 1 = 27 configurações de novo
leaves_search = GridSearchCV(model(), leaves_grid, cv=cv, scoring='roc_auc').fit(X_dev, y_dev)
print(f"leaves grid: best ROC AUC {leaves_search.best_score_:.4f} · {leaves_search.best_params_}")
res = pd.DataFrame(leaves_search.cv_results_)
print(res.groupby('param_max_leaf_nodes')['mean_test_score'].max().round(4).to_string())

leaves grid: best ROC AUC 0.9504 · {'learning_rate': 0.3, 'max_leaf_nodes': 96, 'min_samples_leaf': 5}
param_max_leaf_nodes
16     0.9411
24     0.9436
32     0.9467
48     0.9475
64     0.9491
96     0.9504
128    0.9498
192    0.9487
256    0.9502

Passo 5 — busca aleatória

Preveja. A busca aleatória recebe as mesmas 27 avaliações, mas sobre cinco botões, dois dos quais a grade nem testou. Ela vai bater a grade com folga, empatar ou perder?

from sklearn.model_selection import RandomizedSearchCV

random_search = RandomizedSearchCV(model(), space, n_iter=BUDGET, cv=cv, scoring='roc_auc',
                                   random_state=0).fit(X_dev, y_dev)
print(f"random: best ROC AUC {random_search.best_score_:.4f} · {show(random_search.best_params_)}")

random: best ROC AUC 0.9499 · l2_regularization=8.21, learning_rate=0.228, max_features=0.75, max_leaf_nodes=58, min_samples_leaf=2

A busca aleatória chega a 0,9499, empatada com a grade, com uma configuração que nenhuma grade feita à mão conteria (folhas de tamanho 2, penalidade L2 de 8,21, 75% dos atributos por divisão). Além disso, explorou cinco botões pelo preço que a grade pagou por três.

Simulação: o vencedor em folds novos. A busca aleatória escolheu o vencedor porque ele marcou 0,9499 nos três folds de cv. Pontue de novo as 5 melhores configurações da busca em outro embaralhamento das mesmas linhas, StratifiedKFold(n_splits=3, shuffle=True, random_state=0). O vencedor continua em primeiro? O escore dele cai mais que o das outras?

fresh = StratifiedKFold(n_splits=3, shuffle=True, random_state=0)    # mesmas linhas, outros folds
top5 = pd.DataFrame(random_search.cv_results_).sort_values('rank_test_score').head(5)
for params, s in zip(top5['params'], top5['mean_test_score']):
    again = cross_val_score(model(**params), X_dev, y_dev, cv=fresh, scoring='roc_auc').mean()
    print(f"search {s:.4f} → fresh folds {again:.4f} ({again - s:+.4f})  {show(params)}")

search 0.9499 → fresh folds 0.9490 (-0.0009)  l2_regularization=8.21, learning_rate=0.228, max_features=0.75, max_leaf_nodes=58, min_samples_leaf=2
search 0.9495 → fresh folds 0.9505 (+0.0010)  l2_regularization=0.13, learning_rate=0.0923, max_features=0.75, max_leaf_nodes=71, min_samples_leaf=6
search 0.9485 → fresh folds 0.9444 (-0.0041)  l2_regularization=0.0631, learning_rate=0.11, max_features=0.75, max_leaf_nodes=43, min_samples_leaf=2
search 0.9476 → fresh folds 0.9477 (+0.0001)  l2_regularization=0.00256, learning_rate=0.0226, max_features=0.5, max_leaf_nodes=37, min_samples_leaf=4
search 0.9475 → fresh folds 0.9447 (-0.0027)  l2_regularization=0.157, learning_rate=0.164, max_features=0.75, max_leaf_nodes=33, min_samples_leaf=8

Para fechar a parte fundamental. Em uma frase cada: (1) por que a busca aleatória explorou cinco botões pelo preço que a grade pagou por três? (2) O que um vencedor no canto de uma grade manda você fazer em seguida? (3) Por que o melhor escore que uma busca relata não é o escore do modelo que ela escolheu?

Sua vez (fundamental)

Em casa:

1. A sorte da semente. Rode de novo o passo 5 com random_state=1, 2, 3, 4. Tabele o melhor escore de cada rodada. A variação entre sementes é maior ou menor que o desvio-padrão entre folds dos valores padrão (0,0042)?
2. Uniforme em vez de log-uniforme. Copie space, troque a taxa de aprendizado por uniform(0.01, 0.49) (de scipy.stats) e sorteie 27 configurações com o ParameterSampler. Quantas têm taxa de aprendizado abaixo de 0,05, contra o espaço log-uniforme? Rode de novo o passo 5 com ele: o melhor escore muda?
3. O dobro do orçamento. Rode de novo o passo 5 com n_iter=54. O melhor escore se move mais que o desvio-padrão entre folds? Usando a seção 5, explique por que um best_score_ maior não provaria, sozinho, um modelo melhor.

Abordagem avançada

A parte fundamental gastou um orçamento com quatro métodos e conferiu o que saiu. Esta parte abre os métodos mais espertos e alarga a busca: quando o substituto da otimização bayesiana erra, como o TPE o substitui e como o Optuna escreve um espaço ramificado como código, como parar cedo uma única tentativa e dividir a aposta do halving sucessivo e o que os sistemas de AutoML que buscam o pipeline inteiro acrescentam (e ainda deixam com você). Os blocos de código continuam os das seções 1 a 5.

Esta parte segue um caminho em quatro etapas, desenhado abaixo: 1. questionar o substituto; 2. modelar as boas tentativas; 3. podar e dividir a aposta; 4. o pipeline inteiro. Cada etapa leva à sua seção e, nos slides, uma faixa embaixo mostra em que etapa você está.

Slides da abordagem avançada

Clique nos slides e use as setas ← → para navegar; F abre em tela cheia e Esc mostra todos os slides. Abrir os slides numa aba nova

7. Quando o substituto erra

Etapa 1. As quatro iterações da seção 3 acharam o pico alto na quinta avaliação. Foi o método ou a sorte de três pontos de partida? E o que o substituto supôs para chegar lá?

As hipóteses do substituto importam

Repita a rodada da seção 3 a partir de 100 trios iniciais aleatórios e compare com a busca aleatória, todas com 30 avaliações:

rng = np.random.default_rng(0)
first_hit = {'BO, length scale 0.15': [], 'BO, length scale 0.10': [], 'random search': []}
for run in range(100):
    start = rng.uniform(0, 1, 3)                                  # 3 avaliações aleatórias para começar
    runs = {'BO, length scale 0.15': bayes_opt(start, 27, 0.15)[1],
            'BO, length scale 0.10': bayes_opt(start, 27, 0.10)[1],
            'random search': score(rng.uniform(0, 1, 30))}
    for name, ys in runs.items():                                 # 30 avaliações cada
        hits = np.flatnonzero(ys >= 0.81)                         # a menos de 0.01 do ótimo
        first_hit[name].append(hits[0] + 1 if len(hits) else np.inf)
for name, n in first_hit.items():
    n = np.array(n)
    print(f"{name:22s} reached 0.81 within 10 / 20 / 30 evaluations: "
          f"{np.mean(n <= 10):.0%} / {np.mean(n <= 20):.0%} / {np.mean(n <= 30):.0%}")

BO, length scale 0.15  reached 0.81 within 10 / 20 / 30 evaluations: 40% / 65% / 65%
BO, length scale 0.10  reached 0.81 within 10 / 20 / 30 evaluations: 61% / 92% / 92%
random search          reached 0.81 within 10 / 20 / 30 evaluations: 37% / 54% / 65%

Com escala de comprimento 0,15, a otimização bayesiana está à frente da busca aleatória com 10 avaliações (40% das rodadas contra 37%)... e depois para de melhorar: 35% das rodadas ainda não encontraram o pico alto depois de 30 avaliações. O substituto tem certeza demais de que o escore é suave. Depois que o pico largo em 0,25 foi bem explorado e alguns pontos à direita voltaram com 0,70, ele prevê um 0,70 plano e certo sobre o pico estreito (largura 0,06) e não põe melhoria esperada nenhuma ali. Com escala 0,10, mais perto da largura real do pico, 92% das rodadas o encontram em até 20 avaliações, contra 54% da busca aleatória.

O substituto, portanto, é um modelo como outro qualquer: um modelo errado toma decisões erradas com confiança. Implementações reais não fixam a escala de comprimento à mão como fizemos: elas ajustam os parâmetros do kernel às avaliações por máxima verossimilhança a cada iteração; muitas acrescentam também um pouco de exploração forçada. Elas também precisam lidar com o que um GP trata mal: hiperparâmetros inteiros, categóricos e condicionais, além de espaços com mais de uma dúzia de dimensões. Foi para isso que o próximo método foi criado.

8. TPE e Optuna

Etapa 2. Um processo gaussiano lida mal com inteiros, categorias, botões condicionais e mais de uma dúzia de dimensões. O método que a maioria das bibliotecas usa no lugar dele inverte o modelo.

TPE: uma densidade para as boas tentativas, outra para o resto

O Tree-structured Parzen Estimator (Bergstra, Bardenet, Bengio & Kégl, 2011) inverte a pergunta. Em vez de modelar o escore em função de \(x\), ele modela \(x\) dado o escore. Separe as tentativas até agora nas boas (a melhor fração \(\gamma\)) e no resto; estime uma densidade \(\ell(x)\) dos valores bons e uma densidade \(g(x)\) dos outros, cada uma como a média de uma lombada gaussiana por tentativa (um estimador de Parzen, a estimativa de densidade por kernel da EDA). Depois, amostre onde \(\ell(x)/g(x)\) é maior: valores parecidos com os das tentativas boas e diferentes dos das ruins. Bergstra et al. mostram que maximizar essa razão é o mesmo que maximizar a melhoria esperada sob esse modelo.

Com 12 tentativas aleatórias do mesmo escore, tomando o melhor quarto como "bom" (3 tentativas):

rng = np.random.default_rng(6)
x_tried = rng.uniform(0, 1, 12)                              # 12 tentativas aleatórias até agora
y_tried = score(x_tried)
good = y_tried >= np.quantile(y_tried, 0.75)                 # o melhor quarto: 3 tentativas


def parzen(x, centers, h=0.08):
    """Estimativa de densidade: a média de uma lombada gaussiana por tentativa."""
    return norm.pdf((np.reshape(x, (-1, 1)) - centers) / h).mean(axis=1) / h


l, g = parzen(xs, x_tried[good]), parzen(xs, x_tried[~good])
print('good trials:', ', '.join(f"x = {x:.3f} (score {v:.3f})" for x, v in zip(x_tried[good], y_tried[good])))
for x in [0.1, 0.25, 0.5, 0.7, 0.9]:
    print(f"x = {x}: l = {parzen(x, x_tried[good])[0]:.2f}, g = {parzen(x, x_tried[~good])[0]:.2f}, "
          f"l / g = {parzen(x, x_tried[good])[0] / parzen(x, x_tried[~good])[0]:.2f}")
print(f"next trial: x = {xs[np.argmax(l / g)]:.3f} (score {score(xs[np.argmax(l / g)]):.3f})")

good trials: x = 0.674 (score 0.767), x = 0.330 (score 0.742), x = 0.680 (score 0.777)
x = 0.1: l = 0.03, g = 1.00, l / g = 0.03
x = 0.25: l = 1.01, g = 0.81, l / g = 1.24
x = 0.5: l = 0.46, g = 1.02, l / g = 0.45
x = 0.7: l = 3.19, g = 0.56, l / g = 5.73
x = 0.9: l = 0.07, g = 0.76, l / g = 0.09
next trial: x = 0.710 (score 0.817)

Duas das três tentativas boas estão perto de 0,68, então \(\ell\) é alta ali (3,19) enquanto as ruins são esparsas (\(g = 0{,}56\)): a razão é 5,73 e a próxima tentativa vai para 0,710, que marca 0,817. Perto de 0,1, as tentativas ruins dominam e a razão é 0,03. O TPE segue a densidade das tentativas boas: com outro conjunto de 12 tentativas aleatórias em que as boas se agrupassem em volta do pico largo em 0,25, ele iria para lá primeiro (experimente default_rng(3)).

O que torna o TPE prático é a "árvore" do nome. Cada hiperparâmetro ganha o seu próprio par de densidades unidimensionais, então inteiros, categorias e escalas log são fáceis; um hiperparâmetro que só existe em algumas tentativas é modelado só a partir delas. Não há matriz para inverter, então mil tentativas não custam nada.

Optuna: o espaço de busca é código

O Optuna (Akiba et al., 2019) é a biblioteca que a maioria usa para isso, com o TPE como amostrador padrão (as 10 primeiras tentativas dele são aleatórias, para semear as densidades; o conjunto "bom" é o melhor 10%). Você escreve uma função que recebe um trial, pede valores a ele e devolve o escore. Primeiro, um problema cuja resposta conhecemos: seis hiperparâmetros, dos quais um importa muito, um importa pouco e quatro não importam nada, como na seção grade contra aleatória. O melhor escore possível é 0,870.

import optuna
optuna.logging.set_verbosity(optuna.logging.WARNING)       # silencia o log de cada tentativa


def objective(trial):
    """Seis hiperparâmetros: x0 importa muito, x1 pouco, x2 a x5 nada."""
    x = [trial.suggest_float(f'x{j}', 0, 1) for j in range(6)]
    return 0.70 + 0.12 * np.exp(-((x[0] - 0.72) / 0.08) ** 2) + 0.05 * np.exp(-((x[1] - 0.3) / 0.2) ** 2)


for name, sampler in [('random', optuna.samplers.RandomSampler), ('TPE', optuna.samplers.TPESampler)]:
    best = []
    for seed in range(20):
        study = optuna.create_study(direction='maximize', sampler=sampler(seed=seed))
        study.optimize(objective, n_trials=60)
        best.append(np.maximum.accumulate([t.value for t in study.trials]))
    best = np.array(best)
    print(f"{name:6s} mean best after 10 / 20 / 40 / 60 trials: "
          + ' / '.join(f"{best[:, n - 1].mean():.3f}" for n in (10, 20, 40, 60)))

random mean best after 10 / 20 / 40 / 60 trials: 0.812 / 0.819 / 0.830 / 0.843
TPE    mean best after 10 / 20 / 40 / 60 trials: 0.812 / 0.847 / 0.857 / 0.860

Depois de 10 tentativas os dois são idênticos (0,812): as 10 primeiras tentativas do TPE são aleatórias. A partir daí o TPE se adianta: 0,847 contra 0,819 depois de 20 tentativas; depois de 40, ele está num nível (0,857) que a busca aleatória não alcançou depois de 60 (0,843). O ganho vem de aprender, tentativa após tentativa, qual região de \(x_0\) é boa.

Num problema real, a definição por código compensa de uma segunda maneira: o espaço pode se ramificar. Aqui o Optuna escolhe entre duas famílias para o kc1, cada uma com os próprios botões. Para a regressão logística, as métricas de código são contagens com caudas longas à direita, então passam por log e padronização primeiro:

from sklearn.base import clone
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)


def build(trial):
    """Define-by-run: os hiperparâmetros que existem dependem da família sorteada."""
    family = trial.suggest_categorical('family', ['logistic', 'boosting'])
    if family == 'logistic':
        C = trial.suggest_float('C', 1e-3, 1e2, log=True)
        return make_pipeline(FunctionTransformer(np.log1p), StandardScaler(),
                             LogisticRegression(C=C, max_iter=2000))
    return HistGradientBoostingClassifier(
        learning_rate=trial.suggest_float('learning_rate', 0.01, 0.5, log=True),
        max_leaf_nodes=trial.suggest_int('max_leaf_nodes', 4, 64, log=True),
        min_samples_leaf=trial.suggest_int('min_samples_leaf', 5, 200, log=True),
        random_state=0)


def objective(trial):
    model, scores = build(trial), []
    for fold, (tr, va) in enumerate(cv.split(X_dev, y_dev)):
        m = clone(model).fit(X_dev.iloc[tr], y_dev.iloc[tr])
        scores.append(roc_auc_score(y_dev.iloc[va], m.predict_proba(X_dev.iloc[va])[:, 1]))
        trial.report(np.mean(scores), fold)         # a média corrente depois de cada fold
        if trial.should_prune():                    # pior que a tentativa mediana neste fold?
            raise optuna.TrialPruned()
    return np.mean(scores)


study = optuna.create_study(direction='maximize', sampler=optuna.samplers.TPESampler(seed=0),
                            pruner=optuna.pruners.MedianPruner(n_startup_trials=5, n_warmup_steps=1))
study.optimize(objective, n_trials=40)
trials = study.trials_dataframe()
print(trials.groupby('params_family')['state'].value_counts().unstack(fill_value=0))
print(f"folds fitted: {sum(len(t.intermediate_values) for t in study.trials)} of {40 * 5}")
print(f"best: {study.best_value:.3f} · " + ', '.join(f"{k}={v:.3g}" if isinstance(v, float) else f"{k}={v}"
                                             for k, v in study.best_params.items()))

state          COMPLETE  PRUNED
params_family
boosting              5       4
logistic             25       6
folds fitted: 171 of 200
best: 0.798 · family=logistic, C=0.0297

build é uma função Python comum com um if: quando a tentativa sorteia 'logistic', os botões do boosting nunca são pedidos e o TPE modela C só a partir das tentativas que o tiveram. É isso que "define-by-run" quer dizer. É assim que espaços condicionais, impossíveis de escrever como uma única grade, viram uma linha cada. O resultado também merece leitura. A melhor configuração é uma regressão logística com C ≈ 0,03, com 0,798; o boosting escolhido à mão na primeira seção marcou de 0,795 a 0,806 conforme o embaralhamento, então as duas famílias estão empatadas. O TPE percebeu cedo que as tentativas logísticas eram consistentemente boas e gastou 31 das 40 tentativas nelas. Um espaço de busca que deixa de fora o modelo simples nunca vai poder lhe dizer que você não precisava do complicado.

9. Podar e dividir a aposta

Etapa 3. O halving sucessivo (seção 4) corta cedo as configurações ruins dando primeiro a todas um orçamento pequeno. A mesma aposta pode ser feita uma tentativa de cada vez, dentro de um estudo do Optuna; o risco dela pode ser dividido entre vários torneios.

Poda: pare cedo uma tentativa ruim

O objetivo do kc1 na seção 8 informa a média corrente depois de cada fold (trial.report) e pergunta trial.should_prune(). O podador pela mediana (MedianPruner) responde sim quando o escore corrente da tentativa é pior do que a mediana das tentativas anteriores no mesmo fold; a tentativa é então abandonada. n_startup_trials=5 espera cinco tentativas completas antes de podar qualquer coisa; n_warmup_steps=1 nunca julga uma tentativa só pelo primeiro fold. Aqui 10 das 40 tentativas foram podadas e 171 folds foram ajustados em vez de 200: 15% do orçamento economizado. A economia cresce com o custo dos passos: uma rede neural que informa a perda de validação a cada época pode ser parada depois de 3 épocas de 100.

Podar é apostar que uma tentativa que começa mal termina mal: a aposta do halving sucessivo na seção 4, feita uma tentativa de cada vez. A seção 4 também mediu quando essa aposta erra; o Hyperband a divide.

Hyperband: vários torneios

Quão agressivo deve ser o torneio? Começar 27 configurações com 10 árvores explora amplamente, mas confia num orçamento muito pequeno; começar 4 configurações com o orçamento completo não confia em nada e explora pouco. O Hyperband (Li et al., 2018) não escolhe: ele roda vários torneios (brackets) de agressividade diferente, cada um custando mais ou menos o mesmo. Com um orçamento completo de \(R = 27\) unidades e \(\eta = 3\):

R, eta = 27, 3                          # o orçamento completo, em unidades do menor; o fator
s_max = round(np.log(R) / np.log(eta))  # 3
for s in range(s_max, -1, -1):
    n = int(np.ceil((s_max + 1) / (s + 1) * eta ** s))     # configurações com que o bracket começa
    r = R // eta ** s                                       # orçamento de cada uma na rodada 1
    rounds = [(n // eta ** i, r * eta ** i) for i in range(s + 1)]
    print(f"bracket s = {s}:  " + '  →  '.join(f"{ni:2d} × {ri:2d}" for ni, ri in rounds)
          + f"   cost {sum(ni * ri for ni, ri in rounds)}")

bracket s = 3:  27 ×  1  →   9 ×  3  →   3 ×  9  →   1 × 27   cost 108
bracket s = 2:  12 ×  3  →   4 ×  9  →   1 × 27   cost 99
bracket s = 1:   6 ×  9  →   2 × 27   cost 108
bracket s = 0:   4 × 27   cost 108

O bracket 3 é o torneio agressivo da seção 4 (27 configurações, 1 unidade cada). O bracket 0 é busca aleatória pura: 4 configurações com o orçamento completo, sem eliminação. No meio, os brackets 2 e 1 dividem a aposta. Cada bracket custa cerca de 100 unidades; o conjunto, 423 unidades (umas 16 avaliações completas) para 49 configurações. Se a ordem inicial for confiável, o bracket 3 acha o vencedor barato; se for enganosa, como com árvores e taxas de aprendizado na seção 4, o bracket 0 ainda dá uma chance justa a quem aprende devagar. O BOHB (Falkner, Klein & Hutter, 2018) troca os sorteios aleatórios do Hyperband por densidades à maneira do TPE; o HyperbandPruner do Optuna leva a mesma ideia para a poda.

10. AutoML de pipeline completo

Etapa 4. Tudo até aqui ajustou os botões de uma ou duas famílias de modelos. A forma mais ampla de AutoML busca o pipeline inteiro de uma vez: imputação, codificação, escala, seleção de atributos, a família do modelo, os hiperparâmetros dela e, por fim, como combinar os melhores modelos. Thornton et al. (2013) chamaram isso de problema CASH (Combined Algorithm Selection and Hyperparameter optimization): é o espaço condicional do exemplo do Optuna, com muito mais ramos. Os principais sistemas diferem no que buscam e em como buscam:

Duas ideias se repetem. Ensembles: como uma busca avalia centenas de modelos, dá para fazer a média das previsões de vários modelos bons e diferentes, o que costuma ser melhor que o melhor sozinho; o ensemble selection acrescenta, de forma gulosa e com reposição, o modelo que mais melhora o escore de validação da média. Partidas a quente: uma busca que começa por configurações que funcionaram em problemas parecidos desperdiça menos tentativas. O preço é uma ferramenta muito pesada: o AutoGluon instala vários gigabytes de dependências e treina dezenas de modelos; o auto-sklearn só dá suporte oficial ao Linux. As interfaces são curtas; um fragmento ilustrativo, não executado aqui:

# AutoGluon: pip install autogluon  (pesado: vários GB de dependências)
from autogluon.tabular import TabularPredictor
predictor = TabularPredictor(label='defective', eval_metric='roc_auc').fit(
    train_df, time_limit=600, presets='best_quality')     # faz bagging e empilha muitas famílias
predictor.leaderboard(test_df)

# FLAML: pip install flaml
from flaml import AutoML
automl = AutoML()
automl.fit(X_dev, y_dev, task='classification', metric='roc_auc', time_budget=60)
print(automl.best_estimator, automl.best_config)

Repare no que você ainda precisa escrever: o rótulo, a métrica, o limite de tempo e, antes de tudo isso, train_df e test_df. A seção 5 mostrou o que qualquer busca faz com um vazamento nessas tabelas e com um orçamento longo demais; o resto do trabalho vem abaixo.

O resto do trabalho

O que sobra, mesmo com a melhor busca:

• Formular o problema. O que é previsto, para quem, em que momento e quanto custa um erro. A busca não tem como saber que um alarme falso custa uma hora de um engenheiro e um defeito não detectado custa uma missão.
• Escolher a métrica. A busca maximiza a métrica que você nomeia. No kc1, 84,5% dos módulos são limpos, então uma busca por acurácia aprovaria contente um modelo que nunca aponta nada (veja Métricas de Classificação).
• O desenho da validação. Grupos, ordem no tempo e vazamentos, como na seção 5. Uma divisão aleatória de dados ordenados no tempo será explorada pela busca, não corrigida por ela.
• Qualidade dos dados e atributos. Uma hora gasta num bom atributo costuma valer mais que um dia de busca; em Gradient Boosting, remover uma coluna vazada mudou a precisão média em 0,18; o ajuste de hiperparâmetros, em 0,007.
• Implantação e monitoramento. Um ensemble empilhado de 40 modelos pode ser o mais acurado e o mais difícil de servir, explicar e retreinar. Veja MLOps.

Lab avançado — quatro buscas, um orçamento

O lab fundamental terminou num empate: uma grade feita à mão e a busca aleatória, as duas perto de 0,950, contra 0,9468 dos valores padrão. Esta parte acrescenta os dois métodos que fazem mais que sortear: o TPE, que aprende onde olhar, e o halving sucessivo, que olha 180 configurações pelo preço de 27 contado em linhas. Depois, o vencedor é avaliado uma única vez no conjunto de teste; por fim, uma validação cruzada aninhada avalia a própria busca. Continue no mesmo notebook, na mesma sessão: os passos abaixo usam X_dev, y_dev, X_test, y_test, cv, model, space, BUDGET, show, baseline, grid_search e random_search dos passos 1 a 5.

Passo 6 — TPE com Optuna

Preveja. As 10 primeiras tentativas do TPE são aleatórias; as outras, guiadas pelas boas. Você espera o melhor escore dele claramente acima do da busca aleatória ou dentro do ruído entre folds (cerca de 0,004)?

import optuna
optuna.logging.set_verbosity(optuna.logging.WARNING)


def objective(trial):
    params = {'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.5, log=True),
              'max_leaf_nodes': trial.suggest_int('max_leaf_nodes', 4, 128, log=True),
              'min_samples_leaf': trial.suggest_int('min_samples_leaf', 2, 100, log=True),
              'l2_regularization': trial.suggest_float('l2_regularization', 1e-3, 10, log=True),
              'max_features': trial.suggest_categorical('max_features', [0.5, 0.75, 1.0])}
    return cross_val_score(model(**params), X_dev, y_dev, cv=cv, scoring='roc_auc').mean()


study = optuna.create_study(direction='maximize', sampler=optuna.samplers.TPESampler(seed=0))
study.optimize(objective, n_trials=BUDGET)
print(f"TPE:    best ROC AUC {study.best_value:.4f} · {show(study.best_params)}")

TPE:    best ROC AUC 0.9510 · l2_regularization=1.23, learning_rate=0.189, max_features=1, max_leaf_nodes=102, min_samples_leaf=2

O TPE chega a 0,9510 com 102 folhas e folhas de 2 linhas: a direção do canto da grade (árvores grandes, folhas pequenas), encontrada aprendendo, não enumerando. O espaço é escrito de novo nas palavras do Optuna; log=True faz o papel de loguniform e de log_int.

Passo 7 — halving sucessivo, mesmo orçamento em linhas

Preveja. O halving olha 180 configurações por mais ou menos o mesmo número de linhas de treino que a busca aleatória gastou com 27. Ele vai achar um escore de CV melhor? Vai levar o mesmo tempo?

from sklearn.experimental import enable_halving_search_cv  # noqa: F401
from sklearn.model_selection import HalvingRandomSearchCV

halving = HalvingRandomSearchCV(model(), space, n_candidates=180, factor=3, min_resources=160,
                                cv=cv, scoring='roc_auc', random_state=1).fit(X_dev, y_dev)
rows = sum(n * r for n, r in zip(halving.n_candidates_, halving.n_resources_))
print('candidates per round:', halving.n_candidates_)
print('rows per candidate:  ', halving.n_resources_)
print(f"budget in rows: halving {rows:,} · random search {BUDGET * len(X_dev):,}")
print(f"fits:           halving {3 * sum(halving.n_candidates_)} · random search {3 * BUDGET}")
print(f"halving: best ROC AUC {halving.best_score_:.4f} · {show(halving.best_params_)}")

candidates per round: [180, 60, 20, 7]
rows per candidate:   [160, 480, 1440, 4320]
budget in rows: halving 116,640 · random search 116,721
fits:           halving 801 · random search 81
halving: best ROC AUC 0.9528 · l2_regularization=2.75, learning_rate=0.0823, max_features=0.5, max_leaf_nodes=114, min_samples_leaf=5

O torneio começa com 180 candidatos de 160 linhas cada (com CV 3-fold, cada ajuste treina em cerca de 107 linhas), mantém 60 para 480 linhas, 20 para 1.440 e 7 para 4 320, quase todas as 4.323 linhas de desenvolvimento. Em linhas, o orçamento é o mesmo da busca aleatória: 116.640 contra 116.721. Em ajustes, não: 801 contra 81. Centenas de ajustes minúsculos carregam, cada um, um custo fixo (discretizar os atributos, preparar a parada antecipada), então linhas iguais não significam tempo igual: na máquina usada para esta página, a grade, a busca aleatória e o TPE levaram uns 20 segundos cada; o halving, uns 45. Por esse preço, o halving olhou 180 configurações em vez de 27 e encontrou o melhor escore de CV dos quatro métodos, 0,9528 (114 folhas, folhas de tamanho 5, taxa de aprendizado 0,082).

Simulação: árvores como orçamento. A seção 4 mostrou, no kc1, que um torneio contado em árvores elimina primeiro quem aprende devagar. Confira no phoneme: rode o mesmo torneio com o número de árvores como orçamento, 27 candidatos de 10 a 270 árvores, sem parada antecipada (ela brigaria com o orçamento) e compare as taxas de aprendizado que chegam à última rodada com o vencedor do passo 7.

trees_model = HistGradientBoostingClassifier(max_iter=270, early_stopping=False, random_state=42)
trees_halving = HalvingRandomSearchCV(trees_model, space, n_candidates=27, factor=3,
                                      resource='max_iter', min_resources=10, max_resources=270,
                                      cv=cv, scoring='roc_auc', random_state=1).fit(X_dev, y_dev)
res = pd.DataFrame(trees_halving.cv_results_)
res['lr'] = [p['learning_rate'] for p in res['params']]
for it, r in res.groupby('iter'):
    print(f"round {it + 1}: {len(r):2d} candidates at {r['n_resources'].iloc[0]:3d} trees, "
          f"learning rates {r['lr'].min():.3f} to {r['lr'].max():.3f}")
print(f"winner: {show(trees_halving.best_params_)}, CV ROC AUC {trees_halving.best_score_:.4f}")
print(f"step 7's winner: CV ROC AUC {halving.best_score_:.4f} (learning rate "
      f"{halving.best_params_['learning_rate']:.3f})")

round 1: 27 candidates at  10 trees, learning rates 0.014 to 0.320
round 2:  9 candidates at  30 trees, learning rates 0.025 to 0.320
round 3:  3 candidates at  90 trees, learning rates 0.082 to 0.115
round 4:  1 candidates at 270 trees, learning rates 0.102 to 0.102
winner: l2_regularization=0.481, learning_rate=0.102, max_features=1, max_iter=270, max_leaf_nodes=110, min_samples_leaf=14, CV ROC AUC 0.9569
step 7's winner: CV ROC AUC 0.9528 (learning rate 0.082)

Passo 8 — o melhor escore até aqui, avaliação por avaliação

Preveja. Acompanhe o melhor escore até aqui, avaliação por avaliação. Qual método você espera que comece pior? Por quê? Qual deve melhorar mais tarde?

import matplotlib.pyplot as plt

runs = {'grid': grid_search.cv_results_['mean_test_score'],
        'random': random_search.cv_results_['mean_test_score'],
        'TPE': [t.value for t in study.trials]}
best_so_far = pd.DataFrame({name: np.maximum.accumulate(v) for name, v in runs.items()},
                           index=pd.RangeIndex(1, BUDGET + 1, name='evaluations'))
print(best_so_far.loc[[1, 5, 10, 15, 20, 27]].round(4))
print(f"halving, same budget in rows: {halving.best_score_:.4f}")

ax = best_so_far.plot(drawstyle='steps-post', figsize=(7, 3.5))
ax.axhline(halving.best_score_, color='gray', ls='--', label='halving (same budget)')
ax.axhline(baseline['boosting, defaults'], color='gray', ls=':', label='defaults')
ax.set_ylabel('best CV ROC AUC so far')
ax.legend()
plt.show()

               grid  random     TPE
evaluations
1            0.9327  0.9475  0.9476
5            0.9464  0.9499  0.9489
10           0.9465  0.9499  0.9489
15           0.9467  0.9499  0.9510
20           0.9496  0.9499  0.9510
27           0.9498  0.9499  0.9510
halving, same budget in rows: 0.9528

(script da figura; ele roda o lab de novo, mais a busca aleatória e o TPE com outras quatro sementes)

Leia a tabela linha a linha. Depois de uma avaliação, a grade está em 0,9327: ela lista as configurações numa ordem fixa e a primeira (taxa de aprendizado 0,03, 8 folhas) é lenta e pequena. A busca aleatória e o TPE começam com um sorteio, em 0,9475. A busca aleatória encontrou a melhor configuração até a quinta avaliação e nunca a superou. As 10 primeiras tentativas do TPE também são aleatórias (0,9489 depois de 10); as tentativas guiadas encontraram 0,9510 antes da décima quinta. A grade lista as células com a taxa de aprendizado mudando mais devagar, então a primeira célula de 128 folhas com taxa 0,1 só veio na avaliação 16; a grade chegou a 0,9496 na avaliação 17. O cv_results_ dela também mostra para onde foi o orçamento: as nove células de 8 folhas marcaram todas entre 0,931 e 0,937, qualquer que fosse a taxa de aprendizado. O botão que mais importava aqui (o número de folhas) recebeu só três valores e um terço da grade foi gasto no pior deles: o argumento de Bergstra & Bengio, em dados reais.

Agora as linhas finas da figura. Com outras quatro sementes, a busca aleatória termina entre 0,9512 e 0,9522 e o TPE entre 0,9519 e 0,9532: a rodada do lab (semente 0) foi a mais azarada das cinco para os dois. Semente a semente, o TPE termina à frente em 4 das 5 rodadas, por 0,001 em média, o que é menos que a variação entre sementes de um mesmo método e um quarto do desvio-padrão entre folds. Os 0,9528 do halving caem dentro da faixa do TPE. Neste problema, então, o espaço (árvores grandes, folhas pequenas) importa mais que o método; todo método ganha de 0,003 a 0,006 sobre os valores padrão, mais ou menos um desvio-padrão entre folds.

Passo 9 — o único teste

Preveja. O vencedor foi escolhido entre 261 configurações pelo escore de CV. No conjunto de teste, ele vai marcar acima, perto ou abaixo do escore de CV? Como distinguir a maldição do vencedor de um conjunto de teste que é só mais difícil?

from sklearn.metrics import roc_auc_score

winners = {'grid': (grid_search.best_score_, grid_search.best_params_),
           'random': (random_search.best_score_, random_search.best_params_),
           'TPE': (study.best_value, study.best_params),
           'halving': (halving.best_score_, halving.best_params_)}
name = max(winners, key=lambda k: winners[k][0])        # a regra: melhor escore de CV, decidida antes
cv_score, params = winners[name]
final = model(**params).fit(X_dev, y_dev)
print(f"chosen: the {name} winner · {show(params)}")
print(f"chosen:   CV ROC AUC {cv_score:.4f} · test ROC AUC "
      f"{roc_auc_score(y_test, final.predict_proba(X_test)[:, 1]):.4f}")

reference = model().fit(X_dev, y_dev)                   # não é candidato: a escolha já foi feita
print(f"defaults: CV ROC AUC {baseline['boosting, defaults']:.4f} · test ROC AUC "
      f"{roc_auc_score(y_test, reference.predict_proba(X_test)[:, 1]):.4f}")

chosen: the halving winner · l2_regularization=2.75, learning_rate=0.0823, max_features=0.5, max_leaf_nodes=114, min_samples_leaf=5
chosen:   CV ROC AUC 0.9528 · test ROC AUC 0.9455
defaults: CV ROC AUC 0.9468 · test ROC AUC 0.9374

A regra foi fixada antes de olhar: a configuração com o melhor escore de CV, qualquer que seja o método que a encontrou. É o vencedor do halving, escolhido entre 27 + 27 + 27 + 180 = 261 configurações. No conjunto de teste ele marca 0,9455, 0,007 abaixo do escore de CV dele. É a maldição do vencedor? Os valores padrão, pontuados no mesmo conjunto de teste só como referência (não eram candidatos e nada é escolhido a partir desta linha), caem mais ou menos o mesmo, de 0,9468 para 0,9374: essas 1.081 linhas de teste são simplesmente um pouco mais difíceis do que os folds de desenvolvimento. Medido contra os valores padrão nas mesmas linhas, o ajuste vale +0,008 no conjunto de teste, em linha com o +0,006 que mostrou na validação cruzada.

Passo 10 — uma estimativa honesta da busca

Preveja. A estimativa aninhada roda a busca aleatória inteira dentro de cada um de 3 folds externos e avalia a escolha dela em dados que a busca nunca viu. Ela vai sair acima ou abaixo do best_score_ da própria busca?

outer = StratifiedKFold(n_splits=3, shuffle=True, random_state=7)
nested = cross_val_score(random_search, X_dev, y_dev, cv=outer, scoring='roc_auc')
print(f"random search, best_score_:     {random_search.best_score_:.4f}")
print(f"random search, nested estimate: {nested.mean():.4f} ± {nested.std():.4f}  {nested.round(4)}")

random search, best_score_:     0.9499
random search, nested estimate: 0.9472 ± 0.0067  [0.9463 0.9394 0.9558]

O cross_val_score aceita o objeto de busca como se fosse um modelo: em cada um dos 3 folds externos, ele roda a busca aleatória inteira (27 configurações × 3 folds internos) na parte de treino externa e avalia a configuração escolhida na parte de validação externa, que a busca nunca viu. É a CV aninhada da aula de validação cruzada, aplicada a uma busca. A estimativa honesta de "busca aleatória com 27 tentativas" é 0,9472 ± 0,0067, abaixo dos 0,9499 que a busca relatou para si mesma: 0,003 de otimismo, do mesmo tamanho que todo o ganho da busca sobre os valores padrão (0,9468). Os folds externos vão de 0,9394 a 0,9558, uma variação muito maior que as diferenças entre os quatro métodos. No phoneme, um modelo ajustado é um pouco melhor que os valores padrão; os únicos números que dizem quanto são os medidos em dados que a busca nunca tocou.

Simulação: os valores padrão nos mesmos folds externos. A estimativa aninhada avaliou a busca aleatória em 3 folds externos. Pontue o model() sem ajuste nos mesmos folds externos e olhe o ganho do ajuste fold a fold. O ganho é maior que a variação dele entre folds?

plain = cross_val_score(model(), X_dev, y_dev, cv=outer, scoring='roc_auc')
print(f"defaults, same outer folds: {plain.mean():.4f} ± {plain.std():.4f}  {plain.round(4)}")
print(f"gain of the search, fold by fold: {(nested - plain).round(4)}  mean {(nested - plain).mean():+.4f}")

defaults, same outer folds: 0.9426 ± 0.0058  [0.9407 0.9368 0.9505]
gain of the search, fold by fold: [0.0057 0.0026 0.0053]  mean +0.0045

Para fechar a aula. Em uma frase cada: (1) por que o halving achou o melhor escore de CV dos quatro métodos e quanto isso custou em tempo? (2) Por que o vencedor perdeu 0,007 no teste, e como os valores padrão mostram que isso não é só a maldição do vencedor? (3) O que a estimativa aninhada mede que o best_score_ não mede?

Sua vez (avançado)

Em casa, com mais tempo:

1. A sorte da semente, para dois métodos. Rode de novo os passos 5 e 6 com random_state=1, 2, 3, 4 e TPESampler(seed=1, 2, 3, 4). Tabele o melhor escore de cada rodada. A diferença entre o TPE e a busca aleatória é maior do que a variação entre sementes de um mesmo método? O que isso diz sobre declarar um vencedor a partir de uma rodada de cada?
2. Poda. Reescreva o objetivo do passo 6 como um laço sobre os folds de cv (com trial.report e trial.should_prune()) e crie o estudo com pruner=optuna.pruners.MedianPruner(n_startup_trials=5, n_warmup_steps=1) (o exemplo do kc1 em Poda mostra como). Rode 60 tentativas. Quantas foram podadas? Quantos folds foram ajustados? O melhor escore mudou?
3. Duas famílias. Acrescente um ramo ao objetivo, family = trial.suggest_categorical('family', ['boosting', 'svm']), em que o SVM é make_pipeline(StandardScaler(), SVC(C=..., gamma=...)) com C e gamma em escala log (o cross_val_score com scoring='roc_auc' usa o decision_function dele). Quantas das 40 tentativas o TPE dá a cada família? Qual vence?

Material de aula

Referências

• Bergstra, J.; Bengio, Y. "Random Search for Hyper-Parameter Optimization." JMLR 13 (2012): 281–305 — busca aleatória contra busca em grade e a baixa dimensionalidade efetiva dos problemas de ajuste. texto completo
• Jones, D. R.; Schonlau, M.; Welch, W. J. "Efficient Global Optimization of Expensive Black-Box Functions." Journal of Global Optimization 13 (1998): 455–492 — um substituto de processo gaussiano e a melhoria esperada. DOI
• Snoek, J.; Larochelle, H.; Adams, R. P. "Practical Bayesian Optimization of Machine Learning Algorithms." NeurIPS 25 (2012): 2951–2959 — otimização bayesiana aplicada ao ajuste de modelos de aprendizado de máquina. arXiv:1206.2944
• Bergstra, J.; Bardenet, R.; Bengio, Y.; Kégl, B. "Algorithms for Hyper-Parameter Optimization." NeurIPS 24 (2011): 2546–2554 — o Tree-structured Parzen Estimator (TPE). texto completo
• Akiba, T.; Sano, S.; Yanase, T.; Ohta, T.; Koyama, M. "Optuna: A Next-generation Hyperparameter Optimization Framework." KDD '19 (2019): 2623–2631 — espaços de busca define-by-run e poda. DOI
• Jamieson, K.; Talwalkar, A. "Non-stochastic Best Arm Identification and Hyperparameter Optimization." AISTATS, PMLR 51 (2016): 240–248 — halving sucessivo para busca de hiperparâmetros. texto completo
• Li, L.; Jamieson, K.; DeSalvo, G.; Rostamizadeh, A.; Talwalkar, A. "Hyperband: A Novel Bandit-Based Approach to Hyperparameter Optimization." JMLR 18, no. 185 (2018): 1–52. texto completo
• Falkner, S.; Klein, A.; Hutter, F. "BOHB: Robust and Efficient Hyperparameter Optimization at Scale." ICML, PMLR 80 (2018): 1437–1446 — Hyperband com um modelo à maneira do TPE. texto completo
• Thornton, C.; Hutter, F.; Hoos, H. H.; Leyton-Brown, K. "Auto-WEKA: Combined Selection and Hyperparameter Optimization of Classification Algorithms." KDD '13 (2013): 847–855 — o problema CASH. DOI
• Feurer, M.; Klein, A.; Eggensperger, K.; Springenberg, J.; Blum, M.; Hutter, F. "Efficient and Robust Automated Machine Learning." NeurIPS 28 (2015): 2962–2970 — o auto-sklearn. texto completo
• Caruana, R.; Niculescu-Mizil, A.; Crew, G.; Ksikes, A. "Ensemble Selection from Libraries of Models." ICML (2004) — o ensemble guloso usado pelo auto-sklearn. DOI
• Olson, R. S.; Bartley, N.; Urbanowicz, R. J.; Moore, J. H. "Evaluation of a Tree-based Pipeline Optimization Tool for Automating Data Science." GECCO '16 (2016): 485–492 — o TPOT. DOI
• Erickson, N.; Mueller, J.; Shirkov, A.; Zhang, H.; Larroy, P.; Li, M.; Smola, A. "AutoGluon-Tabular: Robust and Accurate AutoML for Structured Data." (2020). arXiv:2003.06505
• Wang, C.; Wu, Q.; Weimer, M.; Zhu, E. "FLAML: A Fast and Lightweight AutoML Library." MLSys 3 (2021). arXiv:1911.04706
• Hutter, F.; Kotthoff, L.; Vanschoren, J. (orgs.) Automated Machine Learning: Methods, Systems, Challenges. Springer (2019) — o livro de referência da área, de acesso aberto. DOI
• Sayyad Shirabad, J.; Menzies, T. J. The PROMISE Repository of Software Engineering Databases. University of Ottawa (2005) — os dados kc1 usados nesta aula (OpenML 1067).
• Conjunto phoneme, do projeto de reconhecimento de fala ESPRIT 5516 ROARS, distribuído pelo projeto ELENA — os dados do lab (OpenML 1489).
• Optuna — otimização de hiperparâmetros.

Quiz

O quiz cobre a abordagem fundamental.