¿Cómo motiva o se le ocurren una historia sobre la regresión del proceso gaussiano en un weblog principalmente dedicado al aprendizaje profundo?
Fácil. Como lo demuestra las “guerras” de Twitter aparentemente inevitables y confiables que rodean la IA, nada atrae la atención como la controversia y el antagonismo. Entonces, retrocedamos veinte años y encontremos citas de personas que dicen: “¡Aquí vienen los procesos gaussianos, ¡ya no necesitamos molestarnos con esas redes neuronales difíciles y difíciles de ajustar!” Y hoy, aquí estamos; Todos saben algo sobre el aprendizaje profundo pero ¿quién ha oído hablar de los procesos gaussianos?
Si bien cuentos similares cuentan mucho sobre la historia de la ciencia y el desarrollo de opiniones, preferimos un ángulo diferente aquí. En el prefacio de su libro de 2006 sobre Procesos gaussianos para el aprendizaje automático (Rasmussen y Williams 2005)Rasmussen y Williams dicen que se refieren a las “dos culturas”: las disciplinas de las estadísticas y el aprendizaje automático, respectivamente:
Los modelos de procesos gaussianos en cierto sentido reúnen el trabajo en las dos comunidades.
En esta publicación, ese “en algún sentido” se vuelve muy concreto. Veremos una crimson Keras, definida y entrenada de la manera routine, que tiene una capa de proceso gaussiana para su componente principal. La tarea será una regresión multivariada “easy”.
Como aparte, esta “comunidades reunidas”, o formas de pensar, o estrategias de solución, también lo convierte en una buena caracterización common de la probabilidad de flujo de tensor.
Procesos gaussianos
Un proceso gaussiano es un Distribución sobre funciones, donde los valores de función que muestra son conjuntamente gaussianos – En términos generales, una generalización al infinito del gaussiano multivariante. Además del libro de referencia que ya mencionamos (Rasmussen y Williams 2005)hay una serie de buenas presentaciones en la crimson: ver EG https://distill.pub/2019/visual-exploration-paussian-processes/ o https://peterroelants.github.io/posts/gaussian-process-tutorial/. Y como en todo genial, hay un capítulo sobre procesos gaussianos en el difunto David Mackay’s (Mackay 2002) libro.
En esta publicación, utilizaremos la probabilidad de TensorFlow Proceso gaussiano variacional (VGP) Capa, diseñada para trabajar de manera eficiente con “Massive Information”. Como la regresión del proceso gaussiana (GPR, de ahora en adelante) implica la inversión de una matriz de covarianza, posiblemente grande,, se han hecho intentos para diseñar versiones aproximadas, a menudo basadas en principios variacionales. La implementación de TFP se basa en documentos de Titsias (2009) (Titsias 2009) y Hensman et al. (2013) (Hensman, Fusi y Lawrence 2013). En lugar de (p ( mathbf {y} | mathbf {x}) )la probabilidad actual de los datos objetivo dada la entrada actual, trabajamos con una distribución variacional (q ( mathbf {u}) ) Eso actúa como un límite inferior.
Aquí ( Mathbf {u} ) son los valores de función en un conjunto de llamados inducir puntos de índice especificado por el usuario, elegido para cubrir bien el rango de los datos reales. Este algoritmo es mucho más rápido que el GPR “regular”, ya que solo la matriz de covarianza de ( Mathbf {u} ) tiene que ser invertido. Como veremos a continuación, al menos en este ejemplo (así como en otros no descritos aquí) parece ser bastante robusto en cuanto a la cantidad de Puntos inductores aprobado.
Comencemos.
El conjunto de datos
El Conjunto de datos de resistencia a la compresión de concreto es parte del repositorio de aprendizaje automático de UCI. Su página net cube:
El concreto es el materials más importante en la ingeniería civil. La resistencia a la compresión de concreto es una función altamente no lineal de la edad y los ingredientes.
Función altamente no lineal – ¿No suena intrigante? En cualquier caso, debe constituir un caso de prueba interesante para GPR.
Aquí hay una primera mirada.
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 <dbl> 540.0, 540.0, 332.5, 332.5, 198.6, 266.0, 380.0, 380.0, …
$ blast_furnace_slag <dbl> 0.0, 0.0, 142.5, 142.5, 132.4, 114.0, 95.0, 95.0, 114.0,…
$ fly_ash <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ water <dbl> 162, 162, 228, 228, 192, 228, 228, 228, 228, 228, 192, 1…
$ superplasticizer <dbl> 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 <dbl> 1040.0, 1055.0, 932.0, 932.0, 978.4, 932.0, 932.0, 932.0…
$ fine_aggregate <dbl> 676.0, 676.0, 594.0, 594.0, 825.5, 670.0, 594.0, 594.0, …
$ age <dbl> 28, 28, 270, 365, 360, 90, 365, 28, 28, 28, 90, 28, 270,…
$ energy <dbl> 79.986111, 61.887366, 40.269535, 41.052780, 44.296075, 4…
No es tan grande, solo un poco más de 1000 filas, pero aún así, tendremos espacio para experimentar con diferentes números de Puntos inductores.
Tenemos ocho predictores, todos numéricos. Con exclusión de age
(en días), estos representan masas (en kg) en un metro cúbico de concreto. La variable objetivo, energy
se mide en megapascales.
Obtengamos una descripción rápida de las relaciones mutuas.
Verificar una posible interacción (una que un laico podría pensar fácilmente), ¿actúa la concentración de cemento de manera diferente en la resistencia del concreto dependiendo de cuánta agua hay en la mezcla?
Para anclar nuestra percepción futura de qué tan bien lo hace VGP para este ejemplo, se ajustamos a un modelo lineal easy, así como uno que involucra interacciones bidireccionales.
# scale predictors right here already, so information are the identical for all fashions
concrete[, 1:8] <- scale(concrete[, 1:8])
# train-test cut up
set.seed(777)
cut up <- initial_split(concrete, prop = 0.8)
practice <- coaching(cut up)
check <- testing(cut up)
# easy linear mannequin with no interactions
fit1 <- lm(energy ~ ., information = practice)
fit1 %>% abstract()
Name:
lm(formulation = 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
Name:
lm(formulation = 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
También almacenamos las predicciones en el conjunto de pruebas, para una comparación posterior.
linreg_preds1 <- fit1 %>% predict(check[, 1:8])
linreg_preds2 <- fit2 %>% predict(check[, 1:8])
examine <-
information.body(
y_true = check$energy,
linreg_preds1 = linreg_preds1,
linreg_preds2 = linreg_preds2
)
Sin requerido más preprocesamiento, el tfdatasets La tubería de entrada termina bien y corta:
create_dataset <- perform(df, batch_size, shuffle = TRUE) {
df <- as.matrix(df)
ds <-
tensor_slices_dataset(record(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 potential selection for batch dimension ...
batch_size <- 64
train_ds <- create_dataset(practice, batch_size = batch_size)
test_ds <- create_dataset(check, batch_size = nrow(check), shuffle = FALSE)
Y en la creación de modelos.
El modelo
La definición del modelo también es corta, aunque hay algunas cosas para expandirse. No ejecute esto todavía:
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))
)
Dos argumentos a layer_variational_gaussian_process()
Necesita algo de preparación antes de que podamos ejecutar esto. Primero, como nos cube la documentación, kernel_provider
debe ser
una instancia de capa equipada con una @Property, que produce un
PositiveSemidefiniteKernel
instancia”.
En otras palabras, la capa VGP envuelve otra capa de Keras que, que, sí mismoenvoltura o agrupa el flujo de tensor Variables
que contiene los parámetros del núcleo.
Podemos hacer uso de reticulate
es nuevo PyClass
constructor para cumplir con los requisitos anteriores. Usando PyClass
podemos heredar directamente de un objeto Python, agregando y/o anulando métodos o campos como nos gusta, y sí, incluso crear una pitón propiedad.
bt <- import("builtins")
RBFKernelFn <- reticulate::PyClass(
"KernelFn",
inherit = tensorflow::tf$keras$layers$Layer,
record(
`__init__` = perform(self, ...) {
kwargs <- record(...)
tremendous()$`__init__`(kwargs)
dtype <- kwargs[["dtype"]]
self$`_amplitude` = self$add_variable(initializer = initializer_zeros(),
dtype = dtype,
identify = 'amplitude')
self$`_length_scale` = self$add_variable(initializer = initializer_zeros(),
dtype = dtype,
identify = '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`)
)
)
)
)
)
El kernel de proceso gaussiano utilizado es uno de varios disponibles en tfp.math.psd_kernels
(psd
Standing for Optimistic Semidefinite), y probablemente el que viene a la mente primero al pensar en GPR: el exponencial cuadradoo exponenciado cuadrático. La versión utilizada en TFP, con hiperparámetros amplitud (a) y escala de longitud ( lambda )es
[k(x,x’) = 2 a exp (frac{- 0.5 (x−x’)^2}{lambda^2}) ]
Aquí el parámetro interesante es la escala de longitud ( lambda ). Cuando tenemos varias características, sus escalas de longitud, como se inducen por el algoritmo de aprendizaje, reflejan su importancia: si, para alguna característica, ( lambda ) es grande, las respectivas desviaciones al cuadrado de la media no importan tanto. La escala de longitud inversa se puede usar para Determinación de relevancia automática (Neal 1996).
La segunda cosa a cuidar es elegir los puntos de índice iniciales. A partir de los experimentos, las opciones exactas no importan tanto, siempre y cuando los datos estén cubiertos con sensatez. Por ejemplo, una forma alternativa que probamos period construir una distribución empírica (tfd_empírico) de los datos, y luego muestra a partir de ellos. Aquí, en su lugar, solo usamos un – innecesario, ciertamente, dada la disponibilidad de pattern
En r – forma elegante de seleccionar observaciones aleatorias de los datos de entrenamiento:
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]
Un punto interesante a tener en cuenta antes de comenzar la capacitación: el cálculo de los parámetros predictivos posteriores implica una descomposición de Cholesky, que podría fallar si, debido a problemas numéricos, la matriz de covarianza ya no es definitiva. Una acción suficiente para tomar en nuestro caso es hacer todos los cálculos utilizando tf$float64
:
Ahora definimos (de verdad, esta vez) y ejecutamos el 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 1 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, un mayor número de Puntos inductores (Probamos 100 y 200) no tuvieron mucho impacto en el rendimiento de la regresión. Tampoco la elección exacta de las constantes de multiplicación (0.1
y 2
) aplicado al núcleo entrenado Variables
(_amplitude
y _length_scale
)
Haga una gran diferencia en el resultado closing.
Predicciones
Generamos predicciones en el conjunto de pruebas y las agregamos al information.body
que contiene las predicciones de los modelos lineales. Al igual que con otras capas de salida probabilística, “las predicciones” son de hecho distribuciones; Para obtener tensores reales, la muestreamos. Aquí, promedimos más de 10 muestras:
Trazamos las predicciones VGP promedio contra la verdad del suelo, junto con las predicciones del modelo lineal easy (cian) y el modelo que incluye interacciones bidireccionales (violeta):
ggplot(examine, 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(examine$y_true), max(examine$y_true)), ylim = c(min(examine$y_true), max(examine$y_true))) +
ylab("predictions") +
theme(side.ratio = 1)

Figura 1: Predicciones versus verdad terrestre para la regresión lineal (sin interacciones; cian), regresión lineal con interacciones de 2 vías (violeta) y VGP (negro).
Además, comparando MSE para los tres conjuntos de predicciones, vemos
Entonces, el VGP supera a ambas líneas de base. Algo más que podríamos estar interesados: ¿Cómo varían sus predicciones? No tanto como quisiéramos, si construimos las estimaciones de incertidumbre solo de ellas. Aquí trazamos las 10 muestras que dibujamos antes:
samples_df <-
information.body(cbind(examine$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(side.ratio = 1)

Figura 2: Predicciones de 10 muestras consecutivas de la distribución VGP.
Discusión: Relevancia de características
Como se mencionó anteriormente, la escala de longitud inversa se puede usar como un indicador de importancia de la característica. Al usar el ExponentiatedQuadratic
solo núcleos, solo habrá una escala de longitud única; En nuestro ejemplo, la inicial dense
Toma de capa de escala (y además, recombinando) las características.
Alternativamente, podríamos envolver el ExponentiatedQuadratic
en FeatureScaled
núcleo.
FeatureScaled
tiene un adicional scale_diag
Parámetro relacionado con exactamente eso: escala de características. Experimentos con FeatureScaled
(y inicial dense
La capa eliminada, para ser “justa”) mostró un rendimiento ligeramente peor, y los aprendidos scale_diag
Los valores variaron bastante de la ejecución a la ejecución. Por esa razón, elegimos presentar el otro enfoque; Sin embargo, incluimos el código para una envoltura FeatureScaled
En caso de que los lectores deseen experimentar con esto:
ScaledRBFKernelFn <- reticulate::PyClass(
"KernelFn",
inherit = tensorflow::tf$keras$layers$Layer,
record(
`__init__` = perform(self, ...) {
kwargs <- record(...)
tremendous()$`__init__`(kwargs)
dtype <- kwargs[["dtype"]]
self$`_amplitude` = self$add_variable(initializer = initializer_zeros(),
dtype = dtype,
identify = 'amplitude')
self$`_length_scale` = self$add_variable(initializer = initializer_zeros(),
dtype = dtype,
identify = 'length_scale')
self$`_scale_diag` = self$add_variable(
initializer = initializer_ones(),
dtype = dtype,
form = 8L,
identify = '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, si todo lo que le importaba period el rendimiento de la predicción, podría usar FeatureScaled
y mantener la inicial dense
Capa de todos modos. Pero en ese caso, probablemente usarías una crimson neuronal, no un proceso gaussiano, de todos modos …
¡Gracias por leer!
Breiman, Leo. 2001. “Modelado estadístico: las dos culturas (con comentarios y una réplica del autor)”. Estadístico. Sci. 16 (3): 199–231. https://doi.org/10.1214/ss/1009213726.
Hensman, James, Nicolo Fusi y Neil D. Lawrence. 2013. “Procesos gaussianos para Massive Information”. Corrección ABS/1309.6835. http://arxiv.org/abs/1309.6835.
Mackay, David JC 2002. Teoría de la información, algoritmos de inferencia y aprendizaje. Nueva York, NY, EE. UU.: Cambridge College Press.
Neal, Radford M. 1996. Aprendizaje bayesiano para redes neuronales. Berlín, Heidelberg: Springer-Verlag.
Rasmussen, Carl Edward y Christopher Ki Williams. 2005. Procesos gaussianos para el aprendizaje automático (cálculo adaptativo y aprendizaje automático). El MIT prensa.
Titsias, Michalis. 2009. “Aprendizaje variacional de inducir variables en procesos gaussianos escasos”. En Actas de la Docubra Conferencia Internacional sobre Inteligencia Synthetic y Estadísticaseditado por David Van Dyk y Max Welling, 5: 567–74. Actas de la investigación del aprendizaje automático. Hilton Clearwater Seashore Resort, Clearwater Seashore, Florida USA: PMLR. http://proceedings.mlr.press/v5/titsias09a.html.