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.

A aula conta isso em cinco etapas, desenhadas abaixo: 1. enquadrar a busca, uma caixa-preta cara e ruidosa e o espaço em que se busca; 2. grade ou aleatória, as duas buscas que não lembram nada; 3. uma busca que aprende, a otimização bayesiana e o TPE, que usam as tentativas passadas para escolher a próxima; 4. cortar tentativas ruins, poda, halving sucessivo e Hyperband; 5. automatizar e conferir, o pipeline inteiro e o que nenhuma busca conserta. Cada etapa leva à sua seção; nos slides, uma faixa embaixo mostra em que etapa você está.

Slides

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

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 mais abaixo escreve essa árvore como um if.

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 de uma seção mais abaixo: 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.

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.

As hipóteses do substituto importam

Aquela rodada foi boa. Repita-a 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.

TPE e Optuna

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.

Poda: pare cedo uma tentativa ruim

Etapa 4. O objetivo acima 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 próxima seção faz essa aposta de forma sistemática e mede quando ela erra.

Halving sucessivo e Hyperband

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)
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:

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, aqui) ou usar parada antecipada e deixar a taxa de aprendizado e o número de árvores serem ajustados juntos.

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 tabela acima (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 acima, 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.

AutoML de pipeline completo

Etapa 5. 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í que começa a próxima seção.

O que o AutoML não consegue fazer

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.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.

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 acima. 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: quatro buscas, um orçamento

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. Quatro buscas recebem o mesmo orçamento, 27 avaliações de uma validação cruzada 3-fold cada: uma grade, busca aleatória, TPE com Optuna e halving sucessivo com o mesmo número de linhas de treino. Depois, o vencedor é avaliado uma única vez no conjunto de teste; por fim, uma validação cruzada aninhada avalia a própria busca.

O lab precisa do Optuna, que não vem instalado no Colab: rode !pip install optuna na primeira célula. Rode os blocos em ordem, numa mesma sessão; toda saída abaixo é o que o código imprime (tudo junto leva de dois a três minutos). Para trabalhar direto nele, baixe o notebook do lab e abra no Jupyter ou no Colab (Arquivo → Fazer upload de notebook). Com outra versão do Optuna, os números do TPE no passo 6 podem mudar um pouco.

Passo 1 — carregar as vogais

# 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

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

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

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.

Passo 5 — busca aleatória

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.

Passo 6 — TPE com Optuna

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

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).

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

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

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

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.

Sua vez

1. A sorte da semente. 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?
4. A unidade de orçamento errada. Rode de novo o passo 7 com o número de árvores como orçamento: model() com early_stopping=False e resource='max_iter', min_resources=10, max_resources=270, n_candidates=27. Quais taxas de aprendizado sobrevivem até a última rodada? Pontue o vencedor com 270 árvores e compare-o com o vencedor do passo 7, como em O que o torneio supõe.

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