Regressão do processo gaussiano com a probabilidade

Como você motiva ou cria uma história sobre a regressão do processo gaussiano em um weblog dedicado principalmente ao aprendizado profundo?

Fácil. Como demonstrado por “guerras” aparentemente inevitáveis ​​e recorrentes do Twitter em torno da IA, nada atrai atenção como controvérsia e antagonismo. Então, vamos voltar vinte anos e encontrar citações de pessoas dizendo: “Aqui vêm os processos gaussianos, não precisamos mais nos preocupar com aquelas redes neurais complicadas e difíceis de ajustar!” E hoje, aqui estamos nós; todo mundo sabe algo sobre aprendizado profundo, mas quem já ouviu falar de processos gaussianos?

Embora contos semelhantes falem muito sobre a história da ciência e o desenvolvimento de opiniões, preferimos um ângulo diferente aqui. No prefácio de seu livro de 2006 sobre Processos gaussianos para aprendizado de máquina (Rasmussen e Williams 2005)Diz Rasmussen e Williams, referindo -se às “duas culturas” – as disciplinas de estatística e aprendizado de máquina, respectivamente:

Os modelos de processos gaussianos, em certo sentido, reúnem trabalho nas duas comunidades.

Neste publish, esse “em certo sentido” fica muito concreto. Veremos uma rede de Keras, definida e treinada da maneira ordinary, que possui uma camada de processo gaussiana para seu principal constituinte. A tarefa será a regressão multivariada “simples”.

Como um aparte, isso “reunindo comunidades” – ou formas de pensar ou estratégias de solução – também é uma boa caracterização geral da probabilidade de tensorflow.

Processos gaussianos

Um processo gaussiano é um Distribuição sobre funções, onde os valores da função que você amostra é gaussiana em conjunto – Grosso falando, uma generalização para a infinidade do gaussiano multivariado. Além do livro de referência que já mencionamos (Rasmussen e Williams 2005)há várias boas apresentações na rede: veja por exemplo https://distill.pub/2019/visual-exploration-gaussian-processes/ ou https://peterroelants.github.io/postss/gaussian-process-tutorial/. E como em tudo authorized, há um capítulo sobre processos gaussianos no falecido David Mackay’s (Mackay 2002) livro.

Nesta postagem, usaremos a probabilidade do tensorflow Processo gaussiano variacional Camada (VGP), projetada para trabalhar com eficiência com “huge information”. Como a regressão do processo gaussiana (GPR, a partir de agora) envolve a inversão de uma matriz de covariância possivelmente grande – foram feitas tentativas de projetar versões aproximadas, geralmente com base em princípios variacionais. A implementação do TFP é baseada em artigos de Titsias (2009) (Titsias 2009) e Hensman et al. (2013) (Hensman, Fusi e Lawrence 2013). Em vez de (p ( mathbf {y} | mathbf {x}) )a probabilidade actual dos dados de destino, dada a entrada actual, trabalhamos com uma distribuição variacional (q ( mathbf {u}) ) que atua como um limite inferior.

Aqui ( mathbf {u} ) são os valores da função em um conjunto dos chamados induzindo pontos de índice Especificado pelo usuário, escolhido para cobrir bem o intervalo dos dados reais. Este algoritmo é muito mais rápido que o GPR “regular”, como apenas a matriz de covariância de ( mathbf {u} ) tem que ser invertido. Como veremos abaixo, pelo menos neste exemplo (assim como em outros não descritos aqui), parece ser bastante robusto quanto ao número de induzindo pontos passou.

Vamos começar.

O conjunto de dados

O Conjunto de dados de resistência à compressão de concreto faz parte do repositório de aprendizado de máquina da UCI. Sua página da net diz:

O concreto é o materials mais importante na engenharia civil. A resistência à compressão de concreto é uma função altamente não linear da idade e ingredientes.

Função altamente não linear – Isso não parece intrigante? De qualquer forma, deve constituir um caso de teste interessante para o GPR.

Aqui está uma primeira olhada.

library(tidyverse)
library(GGally)
library(visreg)
library(readxl)
library(rsample)
library(reticulate)
library(tfdatasets)
library(keras)
library(tfprobability)

concrete <- read_xls(
  "Concrete_Data.xls",
  col_names = c(
    "cement",
    "blast_furnace_slag",
    "fly_ash",
    "water",
    "superplasticizer",
    "coarse_aggregate",
    "fine_aggregate",
    "age",
    "energy"
  ),
  skip = 1
)

concrete %>% glimpse()

Observations: 1,030
Variables: 9
$ cement              540.0, 540.0, 332.5, 332.5, 198.6, 266.0, 380.0, 380.0, …
$ blast_furnace_slag  0.0, 0.0, 142.5, 142.5, 132.4, 114.0, 95.0, 95.0, 114.0,…
$ fly_ash             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ water               162, 162, 228, 228, 192, 228, 228, 228, 228, 228, 192, 1…
$ superplasticizer    2.5, 2.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0…
$ coarse_aggregate    1040.0, 1055.0, 932.0, 932.0, 978.4, 932.0, 932.0, 932.0…
$ fine_aggregate      676.0, 676.0, 594.0, 594.0, 825.5, 670.0, 594.0, 594.0, …
$ age                 28, 28, 270, 365, 360, 90, 365, 28, 28, 28, 90, 28, 270,…
$ energy            79.986111, 61.887366, 40.269535, 41.052780, 44.296075, 4…

Não é tão grande – apenas um pouco mais de 1000 linhas -, mas ainda assim, teremos espaço para experimentar diferentes números de induzindo pontos.

Temos oito preditores, todos numéricos. Com exceção de age (em dias), estes representam massas (em kg) em um metro cúbico de concreto. A variável de destino, energyé medido em megapascais.

Vamos obter uma rápida visão geral dos relacionamentos mútuos.

Verificando uma possível interação (uma em que um leigo poderia pensar facilmente), a concentração de cimento age de maneira diferente na força do concreto, dependendo da quantidade de água na mistura?

cement_ <- lower(concrete$cement, 3, labels = c("low", "medium", "excessive"))
match <- lm(energy ~ (.) ^ 2, information = cbind(concrete(, 2:9), cement_))
abstract(match)

visreg(match, "cement_", "water", gg = TRUE) + theme_minimal()

Para ancorar nossa percepção futura de quão bem o VGP é para este exemplo, ajustamos um modelo linear simples, bem como um envolvendo interações bidirecionais.

# scale predictors right here already, so information are the identical for all fashions
concrete(, 1:8) <- scale(concrete(, 1:8))

# train-test break up 
set.seed(777)
break up <- initial_split(concrete, prop = 0.8)
practice <- coaching(break up)
take a look at <- testing(break up)

# easy linear mannequin with no interactions
fit1 <- lm(energy ~ ., information = practice)
fit1 %>% abstract()

Name:
lm(system = energy ~ ., information = practice)

Residuals:
    Min      1Q  Median      3Q     Max 
-30.594  -6.075   0.612   6.694  33.032 

Coefficients:
                   Estimate Std. Error t worth Pr(>|t|)    
(Intercept)         35.6773     0.3596  99.204  < 2e-16 ***
cement              13.0352     0.9702  13.435  < 2e-16 ***
blast_furnace_slag   9.1532     0.9582   9.552  < 2e-16 ***
fly_ash              5.9592     0.8878   6.712 3.58e-11 ***
water               -2.5681     0.9503  -2.702  0.00703 ** 
superplasticizer     1.9660     0.6138   3.203  0.00141 ** 
coarse_aggregate     1.4780     0.8126   1.819  0.06929 .  
fine_aggregate       2.2213     0.9470   2.346  0.01923 *  
age                  7.7032     0.3901  19.748  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual customary error: 10.32 on 816 levels of freedom
A number of R-squared:  0.627, Adjusted R-squared:  0.6234 
F-statistic: 171.5 on 8 and 816 DF,  p-value: < 2.2e-16

# two-way interactions
fit2 <- lm(energy ~ (.) ^ 2, information = practice)
fit2 %>% abstract()

Name:
lm(system = energy ~ (.)^2, information = practice)

Residuals:
     Min       1Q   Median       3Q      Max 
-24.4000  -5.6093  -0.0233   5.7754  27.8489 

Coefficients:
                                    Estimate Std. Error t worth Pr(>|t|)    
(Intercept)                          40.7908     0.8385  48.647  < 2e-16 ***
cement                               13.2352     1.0036  13.188  < 2e-16 ***
blast_furnace_slag                    9.5418     1.0591   9.009  < 2e-16 ***
fly_ash                               6.0550     0.9557   6.336 3.98e-10 ***
water                                -2.0091     0.9771  -2.056 0.040090 *  
superplasticizer                      3.8336     0.8190   4.681 3.37e-06 ***
coarse_aggregate                      0.3019     0.8068   0.374 0.708333    
fine_aggregate                        1.9617     0.9872   1.987 0.047256 *  
age                                  14.3906     0.5557  25.896  < 2e-16 ***
cement:blast_furnace_slag             0.9863     0.5818   1.695 0.090402 .  
cement:fly_ash                        1.6434     0.6088   2.700 0.007093 ** 
cement:water                         -4.2152     0.9532  -4.422 1.11e-05 ***
cement:superplasticizer              -2.1874     1.3094  -1.670 0.095218 .  
cement:coarse_aggregate               0.2472     0.5967   0.414 0.678788    
cement:fine_aggregate                 0.7944     0.5588   1.422 0.155560    
cement:age                            4.6034     1.3811   3.333 0.000899 ***
blast_furnace_slag:fly_ash            2.1216     0.7229   2.935 0.003434 ** 
blast_furnace_slag:water             -2.6362     1.0611  -2.484 0.013184 *  
blast_furnace_slag:superplasticizer  -0.6838     1.2812  -0.534 0.593676    
blast_furnace_slag:coarse_aggregate  -1.0592     0.6416  -1.651 0.099154 .  
blast_furnace_slag:fine_aggregate     2.0579     0.5538   3.716 0.000217 ***
blast_furnace_slag:age                4.7563     1.1148   4.266 2.23e-05 ***
fly_ash:water                        -2.7131     0.9858  -2.752 0.006054 ** 
fly_ash:superplasticizer             -2.6528     1.2553  -2.113 0.034891 *  
fly_ash:coarse_aggregate              0.3323     0.7004   0.474 0.635305    
fly_ash:fine_aggregate                2.6764     0.7817   3.424 0.000649 ***
fly_ash:age                           7.5851     1.3570   5.589 3.14e-08 ***
water:superplasticizer                1.3686     0.8704   1.572 0.116289    
water:coarse_aggregate               -1.3399     0.5203  -2.575 0.010194 *  
water:fine_aggregate                 -0.7061     0.5184  -1.362 0.173533    
water:age                             0.3207     1.2991   0.247 0.805068    
superplasticizer:coarse_aggregate     1.4526     0.9310   1.560 0.119125    
superplasticizer:fine_aggregate       0.1022     1.1342   0.090 0.928239    
superplasticizer:age                  1.9107     0.9491   2.013 0.044444 *  
coarse_aggregate:fine_aggregate       1.3014     0.4750   2.740 0.006286 ** 
coarse_aggregate:age                  0.7557     0.9342   0.809 0.418815    
fine_aggregate:age                    3.4524     1.2165   2.838 0.004657 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual customary error: 8.327 on 788 levels of freedom
A number of R-squared:  0.7656,    Adjusted R-squared:  0.7549 
F-statistic: 71.48 on 36 and 788 DF,  p-value: < 2.2e-16

Também armazenamos as previsões no conjunto de testes, para comparação posterior.

linreg_preds1 <- fit1 %>% predict(take a look at(, 1:8))
linreg_preds2 <- fit2 %>% predict(take a look at(, 1:8))

evaluate <-
  information.body(
    y_true = take a look at$energy,
    linreg_preds1 = linreg_preds1,
    linreg_preds2 = linreg_preds2
  )

Sem mais pré -processamento mais necessário, o tfdatasets O pipeline de entrada acaba bem e curto:

create_dataset <- perform(df, batch_size, shuffle = TRUE) {
  
  df <- as.matrix(df)
  ds <-
    tensor_slices_dataset(checklist(df(, 1:8), df(, 9, drop = FALSE)))
  if (shuffle)
    ds <- ds %>% dataset_shuffle(buffer_size = nrow(df))
  ds %>%
    dataset_batch(batch_size = batch_size)
  
}

# only one attainable alternative for batch dimension ...
batch_size <- 64
train_ds <- create_dataset(practice, batch_size = batch_size)
test_ds <- create_dataset(take a look at, batch_size = nrow(take a look at), shuffle = FALSE)

E para modelar a criação.

O modelo

A definição do modelo também é curta, embora haja algumas coisas para expandir. Não execute isso ainda:

mannequin <- keras_model_sequential() %>%
  layer_dense(items = 8,
              input_shape = 8,
              use_bias = FALSE) %>%
  layer_variational_gaussian_process(
    # variety of inducing factors
    num_inducing_points = num_inducing_points,
    # kernel for use by the wrapped Gaussian Course of distribution
    kernel_provider = RBFKernelFn(),
    # output form 
    event_shape = 1, 
    # preliminary values for the inducing factors
    inducing_index_points_initializer = initializer_constant(as.matrix(sampled_points)),
    unconstrained_observation_noise_variance_initializer =
      initializer_constant(array(0.1))
  )

Dois argumentos para layer_variational_gaussian_process() Precisa de alguma preparação antes que possamos realmente executar isso. Primeiro, como a documentação nos diz, kernel_provider deve ser

uma instância de camada equipada com um @property, que produz um PositiveSemidefiniteKernel exemplo”.

Em outras palavras, a camada VGP envolve outra camada de Keras que, em sienvolver -se ou pacotes do tensorflow Variables contendo os parâmetros do kernel.

Podemos usar reticulateé novo PyClass construtor para atender aos requisitos acima. Usando PyClasspodemos herdar diretamente de um objeto Python, adicionando e/ou métodos ou campos de substituição como gostamos – e sim, até criar um python propriedade.

bt <- import("builtins")
RBFKernelFn <- reticulate::PyClass(
  "KernelFn",
  inherit = tensorflow::tf$keras$layers$Layer,
  checklist(
    `__init__` = perform(self, ...) {
      kwargs <- checklist(...)
      tremendous()$`__init__`(kwargs)
      dtype <- kwargs(("dtype"))
      self$`_amplitude` = self$add_variable(initializer = initializer_zeros(),
                                            dtype = dtype,
                                            title = 'amplitude')
      self$`_length_scale` = self$add_variable(initializer = initializer_zeros(),
                                               dtype = dtype,
                                               title = 'length_scale')
      NULL
    },
    
    name = perform(self, x, ...) {
      x
    },
    
    kernel = bt$property(
      reticulate::py_func(
        perform(self)
          tfp$math$psd_kernels$ExponentiatedQuadratic(
            amplitude = tf$nn$softplus(array(0.1) * self$`_amplitude`),
            length_scale = tf$nn$softplus(array(2) * self$`_length_scale`)
          )
      )
    )
  )
)

O kernel do processo gaussiano usado é um dos vários disponíveis em tfp.math.psd_kernels (psd defendendo semidefinita positiva), e provavelmente o que vem à mente primeiro ao pensar em GPR: o exponencial ao quadradoou quadrático exponencial. A versão usada no TFP, com hiperparâmetro amplitude (um) e Escala de comprimento ( lambda )é

(okay (x, x ‘) = 2 a exp ( frac {- 0,5 (x-x’)^2} { lambda^2}) )

Aqui, o parâmetro interessante é a escala de comprimento ( lambda ). Quando temos vários recursos, suas escalas de comprimento – conforme induzidas pelo algoritmo de aprendizado – refletem sua importância: se, para algum recurso, ( lambda ) é grande, os respectivos desvios quadrados da média não importam muito. A escala de comprimento inversa pode ser usada para Determinação de relevância automática (Neal 1996).

A segunda coisa a cuidar é escolher os pontos iniciais do índice. A partir de experimentos, as escolhas exatas não importam muito, desde que os dados sejam sensivelmente cobertos. Por exemplo, uma maneira alternativa de tentamos construir uma distribuição empírica (tfd_empirical) a partir dos dados e, em seguida, amostra dele. Aqui, em vez disso, apenas usamos um – desnecessário, reconhecidamente, dada a disponibilidade de pattern de maneira r – extravagante para selecionar observações aleatórias dos dados de treinamento:

num_inducing_points <- 50

sample_dist <- tfd_uniform(low = 1, excessive = nrow(practice) + 1)
sample_ids <- sample_dist %>%
  tfd_sample(num_inducing_points) %>%
  tf$forged(tf$int32) %>%
  as.numeric()
sampled_points <- practice(sample_ids, 1:8)

Um ponto interessante a ser observado antes de começarmos a treinar: o cálculo dos parâmetros preditivos posteriores envolve uma decomposição de Cholesky, que pode falhar se, devido a problemas numéricos, a matriz de covariância não for mais positiva definida. Uma ação suficiente para tomar em nosso caso é fazer todos os cálculos usando tf$float64:

Agora definimos (de verdade, desta vez) e executamos o modelo.

mannequin <- keras_model_sequential() %>%
  layer_dense(items = 8,
              input_shape = 8,
              use_bias = FALSE) %>%
  layer_variational_gaussian_process(
    num_inducing_points = num_inducing_points,
    kernel_provider = RBFKernelFn(),
    event_shape = 1,
    inducing_index_points_initializer = initializer_constant(as.matrix(sampled_points)),
    unconstrained_observation_noise_variance_initializer =
      initializer_constant(array(0.1))
  )

# KL weight sums to at least one for one epoch
kl_weight <- batch_size / nrow(practice)

# loss that implements the VGP algorithm
loss <- perform(y, rv_y)
  rv_y$variational_loss(y, kl_weight = kl_weight)

mannequin %>% compile(optimizer = optimizer_adam(lr = 0.008),
                  loss = loss,
                  metrics = "mse")

historical past <- mannequin %>% match(train_ds,
                         epochs = 100,
                         validation_data = test_ds)

plot(historical past)

Curiosamente, um número maior de induzindo pontos (Tentamos 100 e 200) não tiveram muito impacto no desempenho da regressão. Nem a escolha exata das constantes de multiplicação (0.1 e 2) aplicado ao kernel treinado Variables (_amplitude e _length_scale)

tfp$math$psd_kernels$ExponentiatedQuadratic(
  amplitude = tf$nn$softplus(array(0.1) * self$`_amplitude`),
  length_scale = tf$nn$softplus(array(2) * self$`_length_scale`)
)

Faça muita diferença no resultado remaining.

Previsões

Geramos previsões no conjunto de testes e as adicionamos ao information.body contendo as previsões dos modelos lineares. Como em outras camadas de saída probabilísticas, “as previsões” são de fato distribuições; Para obter tensores reais, provemos deles. Aqui, calculamos a média de mais de 10 amostras:

yhats <- mannequin(tf$convert_to_tensor(as.matrix(take a look at(, 1:8))))

yhat_samples <-  yhats %>%
  tfd_sample(10) %>%
  tf$squeeze() %>%
  tf$transpose()

sample_means <- yhat_samples %>% apply(1, imply)

evaluate <- evaluate %>%
  cbind(vgp_preds = sample_means)

Plotamos as previsões médias de VGP contra a verdade do solo, juntamente com as previsões do modelo linear simples (ciano) e o modelo, incluindo interações bidirecionais (Violet):

ggplot(evaluate, aes(x = y_true)) +
  geom_abline(slope = 1, intercept = 0) +
  geom_point(aes(y = vgp_preds, coloration = "VGP")) +
  geom_point(aes(y = linreg_preds1, coloration = "easy lm"), alpha = 0.4) +
  geom_point(aes(y = linreg_preds2, coloration = "lm w/ interactions"), alpha = 0.4) +
  scale_colour_manual("", 
                      values = c("VGP" = "black", "easy lm" = "cyan", "lm w/ interactions" = "violet")) +
  coord_cartesian(xlim = c(min(evaluate$y_true), max(evaluate$y_true)), ylim = c(min(evaluate$y_true), max(evaluate$y_true))) +
  ylab("predictions") +
  theme(facet.ratio = 1) 

Figura 1: Previsões vs. Verdade do solo para regressão linear (sem interações; ciano), regressão linear com interações bidirecionais (violeta) e VGP (preto).

Além disso, comparando MSEs para os três conjuntos de previsões, vemos

mse <- perform(y_true, y_pred) {
  sum((y_true - y_pred) ^ 2) / size(y_true)
}

lm_mse1 <- mse(evaluate$y_true, evaluate$linreg_preds1) # 117.3111
lm_mse2 <- mse(evaluate$y_true, evaluate$linreg_preds2) # 80.79726
vgp_mse <- mse(evaluate$y_true, evaluate$vgp_preds)     # 58.49689

Portanto, o VGP de fato supera as duas linhas de base. Outra coisa em que podemos estar interessados: como suas previsões variam? Não tanto quanto desejarmos, se formos para construir estimativas de incerteza deles sozinhas. Aqui planejamos as 10 amostras que desenhamos antes:

samples_df <-
  information.body(cbind(evaluate$y_true, as.matrix(yhat_samples))) %>%
  collect(key = run, worth = prediction, -X1) %>% 
  rename(y_true = "X1")

ggplot(samples_df, aes(y_true, prediction)) +
  geom_point(aes(coloration = run),
             alpha = 0.2,
             dimension = 2) +
  geom_abline(slope = 1, intercept = 0) +
  theme(legend.place = "none") +
  ylab("repeated predictions") +
  theme(facet.ratio = 1)

Figura 2: Previsões de 10 amostras consecutivas da distribuição VGP.

Discussão: Relevância do recurso

Como mencionado acima, a escala de comprimento inversa pode ser usada como um indicador de importância do recurso. Ao usar o ExponentiatedQuadratic Kernel sozinho, haverá apenas uma única escala de comprimento; Em nosso exemplo, o inicial dense camada toma escala (e adicionalmente, recombinando) os recursos.

Como alternativa, poderíamos embrulhar o ExponentiatedQuadratic em um FeatureScaled Kernel.
FeatureScaled tem um adicional scale_diag Parâmetro relacionado a exatamente isso: escala de recursos. Experimentos com FeatureScaled (e inicial dense a camada removida, para ser “justa”) mostrou um desempenho um pouco pior, e o aprendido scale_diag Os valores variaram bastante de execução para executar. Por esse motivo, optamos por apresentar a outra abordagem; No entanto, incluímos o código para um embrulho FeatureScaled Caso os leitores gostarem de experimentar com isso:

ScaledRBFKernelFn <- reticulate::PyClass(
  "KernelFn",
  inherit = tensorflow::tf$keras$layers$Layer,
  checklist(
    `__init__` = perform(self, ...) {
      kwargs <- checklist(...)
      tremendous()$`__init__`(kwargs)
      dtype <- kwargs(("dtype"))
      self$`_amplitude` = self$add_variable(initializer = initializer_zeros(),
                                            dtype = dtype,
                                            title = 'amplitude')
      self$`_length_scale` = self$add_variable(initializer = initializer_zeros(),
                                               dtype = dtype,
                                               title = 'length_scale')
      self$`_scale_diag` = self$add_variable(
        initializer = initializer_ones(),
        dtype = dtype,
        form = 8L,
        title = 'scale_diag'
      )
      NULL
    },
    
    name = perform(self, x, ...) {
      x
    },
    
    kernel = bt$property(
      reticulate::py_func(
        perform(self)
          tfp$math$psd_kernels$FeatureScaled(
            kernel = tfp$math$psd_kernels$ExponentiatedQuadratic(
              amplitude = tf$nn$softplus(array(1) * self$`_amplitude`),
              length_scale = tf$nn$softplus(array(2) * self$`_length_scale`)
            ),
            scale_diag = tf$nn$softplus(array(1) * self$`_scale_diag`)
          )
      )
    )
  )
)

Finalmente, se tudo o que você se importava period o desempenho da previsão, você poderia usar FeatureScaled e mantenha o inicial dense camada da mesma forma. Mas, nesse caso, você provavelmente usaria uma rede neural – não um processo gaussiano – de qualquer maneira …

Obrigado pela leitura!