Generalized_Linear_Models.ipynb

Kernel: Python 3

Copyright 2018 The TensorFlow Probability Authors.

Licensed under the Apache License, Version 2.0 (the "License");

#@title Licensed under the Apache License, Version 2.0 (the "License"); { display-mode: "form" }
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

Modelos lineales generalizados

Ver en TensorFlow.org

Ejecutar en Google Colab

Ver fuente en GitHub

Descargar el bloc de notas

En este bloc de notas, se presentan los Modelos Lineales Generalizados (GLM) a través de un ejemplo resuelto. Este ejemplo se resuelve de dos maneras diferentes con dos algoritmos para ajustar los GLM de manera eficiente en TensorFlow Probability: puntuación de Fisher para datos densos y descenso de gradiente proximal en coordenadas para datos dispersos. Se comparan los coeficientes ajustados con los coeficientes verdaderos y, en el caso de un descenso de gradiente proximal en coordenadas, con la salida del algoritmo glmnet similar de R. Finalmente, se ofrecen más detalles matemáticos y derivaciones de varias propiedades clave de los GLM.

Antecedentes

Un modelo lineal generalizado (GLM) es un modelo lineal ($\eta = x^\top \beta$) envuelto en una transformación (función de enlace) y equipado con una distribución de respuesta de una familia exponencial. La elección de la función de enlace y la distribución de la respuesta es muy flexible, lo que otorga una gran expresividad a los GLM. Los detalles completos, incluida una presentación secuencial de todas las definiciones y resultados de los GLM en notación inequívoca, se encuentran en "Derivación de hechos del GLM" a continuación. Resumimos:

En un GLM, una distribución predictiva para la variable de respuesta $Y$ se asocia con un vector de predictores observados $x$. La distribución tiene la siguiente forma:

$$\begin{align*} p(y \, |\, x) &= m(y, \phi) \exp\left(\frac{\theta\, T(y) - A(\theta)}{\phi}\right) \\ \theta &:= h(\eta) \\ \eta &:= x^\top \beta \end{align*}$$

Aquí $\beta$ son los parámetros ("ponderaciones"), $\phi$ un hiperparámetro que representa la dispersión ("varianza") y $m$, $h$, $T$, $A$ se caracterizan por la familia de modelos especificada por el usuario.

La media de $Y$ depende de $x$ por la composición de la respuesta lineal $\eta$ y la función de enlace (inversa), es decir:

$$\mu := g^{-1}(\eta)$$

donde $g$ es la conocida función de enlace. En TFP, la elección de la función de enlace y la familia de modelos se especifican conjuntamente mediante una subclase tfp.glm.ExponentialFamily. Entre los ejemplos, se incluyen:

• tfp.glm.Normal, también conocido como "regresión lineal"

• tfp.glm.Bernoulli, también conocido como "regresión logística"

• tfp.glm.Poisson, también conocido como "regresión de Poisson"

• tfp.glm.BernoulliNormalCDF, también conocido como "regresión probit".

TFP prefiere nombrar familias modelo según la distribución sobre Y en lugar de la función de enlace, ya que las distribuciones tfp.Distribution ya son ciudadanos de primera clase. Si el nombre de la subclase tfp.glm.ExponentialFamily contiene una segunda palabra, esto indica una función de enlace no canónica.

Los GLM tienen varias propiedades notables que permiten la implementación eficiente del estimador de máxima verosimilitud. Entre estas propiedades se destacan las fórmulas simples para el gradiente de la probabilidad logarítmica $\ell$ y para la matriz de información de Fisher, que es el valor esperado del hessiano de la probabilidad logarítmica negativa bajo un nuevo muestreo de la respuesta bajo los mismos predictores. Es decir:

$$\begin{align*} \nabla_\beta\, \ell(\beta\, ;\, \mathbf{x}, \mathbf{y}) &= \mathbf{x}^\top \,\text{diag}\left(\frac{ {\textbf{Mean}_T}'(\mathbf{x} \beta) }{ {\textbf{Var}_T}(\mathbf{x} \beta) }\right) \left(\mathbf{T}(\mathbf{y}) - {\textbf{Mean}_T}(\mathbf{x} \beta)\right) \\ \mathbb{E}_{Y_i \sim \text{GLM} | x_i} \left[ \nabla_\beta^2\, \ell(\beta\, ;\, \mathbf{x}, \mathbf{Y}) \right] &= -\mathbf{x}^\top \,\text{diag}\left( \frac{ \phi\, {\textbf{Mean}_T}'(\mathbf{x} \beta)^2 }{ {\textbf{Var}_T}(\mathbf{x} \beta) }\right)\, \mathbf{x} \end{align*}$$

donde $\mathbf{x}$ es la matriz cuya $i$ésima fila es el vector predictor para la $i$ésima muestra de datos, y $\mathbf{y}$ es el vector cuya $i$ésima coordenada es la respuesta observada para la $i$ésima muestra de datos. Aquí (en términos generales), ${\text{Mean}_T}(\eta) := \mathbb{E}[T(Y),|,\eta]$ y ${\text{Var}_T}(\ eta) := \text{Var}[T(Y),|,\eta]$, y la negrita denota vectorización de estas funciones. Todos los detalles de lo que representan las distribuciones de estas expectativas y varianzas se pueden encontrar en "Derivación de hechos de GLM" a continuación.

Un ejemplo

En esta sección, describimos y mostramos brevemente dos algoritmos de ajuste de GLM integrados en TensorFlow Probability: puntuación de Fisher (tfp.glm.fit) y descenso de gradiente proximal por coordenadas (tfp.glm.fit_sparse).

Conjunto de datos sintéticos

Supongamos que cargamos algún conjunto de datos de entrenamiento.

import numpy as np
import pandas as pd
import scipy
import tensorflow.compat.v2 as tf
tf.enable_v2_behavior()
import tensorflow_probability as tfp
tfd = tfp.distributions

def make_dataset(n, d, link, scale=1., dtype=np.float32):
  model_coefficients = tfd.Uniform(
      low=-1., high=np.array(1, dtype)).sample(d, seed=42)
  radius = np.sqrt(2.)
  model_coefficients *= radius / tf.linalg.norm(model_coefficients)
  mask = tf.random.shuffle(tf.range(d)) < int(0.5 * d)
  model_coefficients = tf.where(
      mask, model_coefficients, np.array(0., dtype))
  model_matrix = tfd.Normal(
      loc=0., scale=np.array(1, dtype)).sample([n, d], seed=43)
  scale = tf.convert_to_tensor(scale, dtype)
  linear_response = tf.linalg.matvec(model_matrix, model_coefficients)
  
  if link == 'linear':
    response = tfd.Normal(loc=linear_response, scale=scale).sample(seed=44)
  elif link == 'probit':
    response = tf.cast(
        tfd.Normal(loc=linear_response, scale=scale).sample(seed=44) > 0,
                   dtype)
  elif link == 'logit':
    response = tfd.Bernoulli(logits=linear_response).sample(seed=44)
  else:
    raise ValueError('unrecognized true link: {}'.format(link))
  return model_matrix, response, model_coefficients, mask

Nota: Conéctese a un entorno de ejecución local.

En este bloc de notas, compartimos datos entre los núcleos de Python y R a través de archivos locales. Para habilitar este uso compartido, utilice tiempos de ejecución en la misma máquina donde tiene permiso para leer y escribir archivos locales.

x, y, model_coefficients_true, _ = [t.numpy() for t in make_dataset(
    n=int(1e5), d=100, link='probit')]

DATA_DIR = '/tmp/glm_example'
tf.io.gfile.makedirs(DATA_DIR)
with tf.io.gfile.GFile('{}/x.csv'.format(DATA_DIR), 'w') as f:
  np.savetxt(f, x, delimiter=',')
with tf.io.gfile.GFile('{}/y.csv'.format(DATA_DIR), 'w') as f:
  np.savetxt(f, y.astype(np.int32) + 1, delimiter=',', fmt='%d')
with tf.io.gfile.GFile(
    '{}/model_coefficients_true.csv'.format(DATA_DIR), 'w') as f:
  np.savetxt(f, model_coefficients_true, delimiter=',')

Sin regularización L1

La función tfp.glm.fit implementa la puntuación de Fisher, que toma como uno de sus argumentos:

• model_matrix = $\mathbf{x}$

• response = $\mathbf{y}$

• model = invocable que, dado el argumento $\boldsymbol{\eta}$, devuelve el triple $\left( {\textbf{Mean}_T}(\boldsymbol{\eta}), {\textbf{Var}_T}(\boldsymbol{\eta}), {\textbf{Mean}_T}'(\boldsymbol{\eta}) \right)$.

Recomendamos que model sea una instancia de la clase tfp.glm.ExponentialFamily. Hay varias implementaciones prediseñadas disponibles, por lo que para los GLM más comunes no se requiere ningún código personalizado.

@tf.function(autograph=False)
def fit_model():
  model_coefficients, linear_response, is_converged, num_iter = tfp.glm.fit(
      model_matrix=x, response=y, model=tfp.glm.BernoulliNormalCDF())
  log_likelihood = tfp.glm.BernoulliNormalCDF().log_prob(y, linear_response)
  return (model_coefficients, linear_response, is_converged, num_iter,
          log_likelihood)
 
[model_coefficients, linear_response, is_converged, num_iter,
 log_likelihood] = [t.numpy() for t in fit_model()]

print(('is_converged: {}\n'
       '    num_iter: {}\n'
       '    accuracy: {}\n'
       '    deviance: {}\n'
       '||w0-w1||_2 / (1+||w0||_2): {}'
      ).format(
    is_converged,
    num_iter,
    np.mean((linear_response > 0.) == y),
    2. * np.mean(log_likelihood),
    np.linalg.norm(model_coefficients_true - model_coefficients, ord=2) /
        (1. + np.linalg.norm(model_coefficients_true, ord=2))
    ))

is_converged: True
    num_iter: 6
    accuracy: 0.75241
    deviance: -0.992436110973
||w0-w1||_2 / (1+||w0||_2): 0.0231555201462

Detalles matemáticos

La puntuación de Fisher es una modificación del método de Newton que nos permite encontrar la estimación de máxima verosimilitud

$$\hat\beta := \underset{\beta}{\text{arg max}}\ \ \ell(\beta\ ;\ \mathbf{x}, \mathbf{y}).$$

El método de Newton básico, que busca ceros en el gradiente del logaritmo de verosimilitud, seguiría la regla de actualización

$$\beta^{(t+1)}_{\text{Newton}} := \beta^{(t)} - \alpha \left( \nabla^2_\beta\, \ell(\beta\ ;\ \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)}}^{-1} \left( \nabla_\beta\, \ell(\beta\ ;\ \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)}}$$

donde $\alpha \in (0, 1]$ es una tasa de aprendizaje que se usa para controlar el tamaño del paso.

En la puntuación de Fisher, el hessiano se reemplaza por la matriz de información negativa de Fisher:

$$\begin{align*} \beta^{(t+1)} &:= \beta^{(t)} - \alpha\, \mathbb{E}_{ Y_i \sim p_{\text{OEF}(m, T)}(\cdot | \theta = h(x_i^\top \beta^{(t)}), \phi) } \left[ \left( \nabla^2_\beta\, \ell(\beta\ ;\ \mathbf{x}, \mathbf{Y}) \right)_{\beta = \beta^{(t)}} \right]^{-1} \left( \nabla_\beta\, \ell(\beta\ ;\ \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)}} $3mm] \end{align*}$$

[Tenga en cuenta que aquí $\mathbf{Y} = (Y_i)_{i=1}^{n}$ es aleatorio, mientras que $\mathbf{y}$ sigue siendo el vector de respuestas observadas].

Mediante las fórmulas en "Ajuste de parámetros del GLM a los datos" a continuación, esto se simplifica de la siguiente manera:

$$\begin{align*} \beta^{(t+1)} &= \beta^{(t)} + \alpha \left( \mathbf{x}^\top \text{diag}\left( \frac{ \phi\, {\textbf{Mean}_T}'(\mathbf{x} \beta^{(t)})^2 }{ {\textbf{Var}_T}(\mathbf{x} \beta^{(t)}) }\right)\, \mathbf{x} \right)^{-1} \left( \mathbf{x}^\top \text{diag}\left(\frac{ {\textbf{Mean}_T}'(\mathbf{x} \beta^{(t)}) }{ {\textbf{Var}_T}(\mathbf{x} \beta^{(t)}) }\right) \left(\mathbf{T}(\mathbf{y}) - {\textbf{Mean}_T}(\mathbf{x} \beta^{(t)})\right) \right). \end{align*}$$

Con regularización L1

tfp.glm.fit_sparse implementa un ajustador del GLM más adecuado para conjuntos de datos dispersos, basado en el algoritmo de Yuan, Ho y Lin 2012. Sus características incluyen lo que sigue:

• Regularización L1

• Sin inversiones de matrices

• Pocas evaluaciones del gradiente y del hessiano

En primer lugar, presentamos un ejemplo de uso del código. Los detalles del algoritmo se detallan mejor en "Detalles del algoritmo para tfp.glm.fit_sparse" a continuación.

model = tfp.glm.Bernoulli()
model_coefficients_start = tf.zeros(x.shape[-1], np.float32)
@tf.function(autograph=False)
def fit_model():
  return tfp.glm.fit_sparse(
    model_matrix=tf.convert_to_tensor(x),
    response=tf.convert_to_tensor(y),
    model=model,
    model_coefficients_start=model_coefficients_start,
    l1_regularizer=800.,
    l2_regularizer=None,
    maximum_iterations=10,
    maximum_full_sweeps_per_iteration=10,
    tolerance=1e-6,
    learning_rate=None)

model_coefficients, is_converged, num_iter = [t.numpy() for t in fit_model()]
coefs_comparison = pd.DataFrame({
  'Learned': model_coefficients,
  'True': model_coefficients_true,
})
  
print(('is_converged: {}\n'
       '    num_iter: {}\n\n'
       'Coefficients:').format(
    is_converged,
    num_iter))
coefs_comparison

is_converged: True
    num_iter: 1

Coefficients:

Tenga en cuenta que los coeficientes aprendidos tienen el mismo patrón de dispersión que los coeficientes verdaderos.

# Save the learned coefficients to a file.
with tf.io.gfile.GFile('{}/model_coefficients_prox.csv'.format(DATA_DIR), 'w') as f:
  np.savetxt(f, model_coefficients, delimiter=',')

Comparación con glmnet de R

Comparamos la salida del descenso de gradiente proximal por coordenadas con la de glmnet de R, que usa un algoritmo similar.

NOTA: Para ejecutar esta sección, se debe cambiar a un tiempo de ejecución de Colab de R.

suppressMessages({
  library('glmnet')
})

data_dir <- '/tmp/glm_example'
x <- as.matrix(read.csv(paste(data_dir, '/x.csv', sep=''),
                        header=FALSE))
y <- as.matrix(read.csv(paste(data_dir, '/y.csv', sep=''),
                        header=FALSE, colClasses='integer'))

fit <- glmnet(
x = x,
y = y,
family = "binomial",  # Logistic regression
alpha = 1,  # corresponds to l1_weight = 1, l2_weight = 0
standardize = FALSE,
intercept = FALSE,
thresh = 1e-30,
type.logistic = "Newton"
)

write.csv(as.matrix(coef(fit, 0.008)),
          paste(data_dir, '/model_coefficients_glmnet.csv', sep=''),
          row.names=FALSE)

Compare R, TFP y coeficientes verdaderos (Nota: Vuelva al núcleo de Python)

DATA_DIR = '/tmp/glm_example'
with tf.io.gfile.GFile('{}/model_coefficients_glmnet.csv'.format(DATA_DIR),
                       'r') as f:
  model_coefficients_glmnet = np.loadtxt(f,
                                   skiprows=2  # Skip column name and intercept
                               )

with tf.io.gfile.GFile('{}/model_coefficients_prox.csv'.format(DATA_DIR),
                       'r') as f:
  model_coefficients_prox = np.loadtxt(f)

with tf.io.gfile.GFile(
    '{}/model_coefficients_true.csv'.format(DATA_DIR), 'r') as f:
  model_coefficients_true = np.loadtxt(f)

coefs_comparison = pd.DataFrame({
    'TFP': model_coefficients_prox,
    'R': model_coefficients_glmnet,
    'True': model_coefficients_true,
})
coefs_comparison

Detalles del algoritmo para tfp.glm.fit_sparse

Presentamos el algoritmo como una secuencia de tres modificaciones al método de Newton. En cada uno, la regla de actualización de $\beta$ se basa en un vector $s$ y una matriz $H$ que se aproximan al gradiente y al hessiano de la probabilidad logarítmica. En el paso $t$, elegimos una coordenada $j^{(t)}$ para cambiar y actualizamos $\beta$ de acuerdo con la regla de actualización:

$$\begin{align*} u^{(t)} &:= \frac{ \left( s^{(t)} \right)_{j^{(t)}} }{ \left( H^{(t)} \right)_{j^{(t)},\, j^{(t)}} } $3mm] \beta^{(t+1)} &:= \beta^{(t)} - \alpha\, u^{(t)} \,\text{onehot}(j^{(t)}) \end{align*}$$

Esta actualización es un paso similar a Newton con una tasa de aprendizaje $\alpha$. Excepto por la pieza final (regularización L1), las siguientes modificaciones difieren solo en cómo actualizan $s$ y $H$.

Punto de partida: método de Newton por coordenadas

En el método de Newton por coordenadas, establecemos $s$ y $H$ en el gradiente verdadero y el hessiano de la probabilidad logarítmica:

$$\begin{align*} s^{(t)}_{\text{vanilla}} &:= \left( \nabla_\beta\, \ell(\beta \,;\, \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)}} \\ H^{(t)}_{\text{vanilla}} &:= \left( \nabla^2_\beta\, \ell(\beta \,;\, \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)}} \end{align*}$$

Menos evaluaciones del gradiente y del hessiano

El gradiente y el hessiano de la probabilidad logarítmica suelen ser costosos de calcular, por lo que a menudo vale la pena aproximarlos. Podemos hacerlo de la siguiente manera:

• Por lo general, se aproxima el hessiano como localmente constante y se aproxima el gradiente al primer orden con ayuda del hessiano (aproximado):

$$\begin{align*} H_{\text{approx}}^{(t+1)} &:= H^{(t)} \\ s_{\text{approx}}^{(t+1)} &:= s^{(t)} + H^{(t)} \left( \beta^{(t+1)} - \beta^{(t)} \right) \end{align*}$$

• De vez en cuando, se ejecuta un paso de actualización "básico" como el anterior, configurando $s^{(t+1)}$ con el gradiente exacto y $H^{(t+1)}$ con el hessiano exacto de la probabilidad logarítmica, evaluado en $\beta^{(t+1)}$.

Sustitución de la información negativa de Fisher por el hessiano

Para reducir aún más el costo de los pasos de actualización básicos, podemos establecer $H$ en la matriz de información negativa de Fisher (que se puede calcular eficientemente si se usan las fórmulas en "Ajuste de parámetros del GLM a los datos" a continuación) en lugar del hessiano exacto:

$$\begin{align*} H_{\text{Fisher}}^{(t+1)} &:= \mathbb{E}_{Y_i \sim p_{\text{OEF}(m, T)}(\cdot | \theta = h(x_i^\top \beta^{(t+1)}), \phi)} \left[ \left( \nabla_\beta^2\, \ell(\beta\, ;\, \mathbf{x}, \mathbf{Y}) \right)_{\beta = \beta^{(t+1)}} \right] \\ &= -\mathbf{x}^\top \,\text{diag}\left( \frac{ \phi\, {\textbf{Mean}_T}'(\mathbf{x} \beta^{(t+1)})^2 }{ {\textbf{Var}_T}(\mathbf{x} \beta^{(t+1)}) }\right)\, \mathbf{x} \\ s_{\text{Fisher}}^{(t+1)} &:= s_{\text{vanilla}}^{(t+1)} \\ &= \left( \mathbf{x}^\top \,\text{diag}\left(\frac{ {\textbf{Mean}_T}'(\mathbf{x} \beta^{(t+1)}) }{ {\textbf{Var}_T}(\mathbf{x} \beta^{(t+1)}) }\right) \left(\mathbf{T}(\mathbf{y}) - {\textbf{Mean}_T}(\mathbf{x} \beta^{(t+1)})\right) \right) \end{align*}$$

Regularización L1 mediante descenso de gradiente proximal

Para incorporar la regularización L1, reemplazamos la regla de actualización.

$$\beta^{(t+1)} := \beta^{(t)} - \alpha\, u^{(t)} \,\text{onehot}(j^{(t)})$$

con la regla de actualización más general

$$\begin{align*} \gamma^{(t)} &:= -\frac{\alpha\, r_{\text{L1}}}{\left(H^{(t)}\right)_{j^{(t)},\, j^{(t)}}} $2mm] \left(\beta_{\text{reg}}^{(t+1)}\right)_j &:= \begin{cases} \beta^{(t+1)}_j &\text{if } j \neq j^{(t)} \\ \text{SoftThreshold} \left( \beta^{(t)}_j - \alpha\, u^{(t)} ,\ \gamma^{(t)} \right) &\text{if } j = j^{(t)} \end{cases} \end{align*}$$

donde ParseError: KaTeX parse error: Expected 'EOF', got '&' at position 15: r_{\text{L1}} &̲gt; 0 es una constante proporcionada (el coeficiente de regularización L1) y $\text{SoftThreshold}$ es el operador de umbral suave, definido por

$$\text{SoftThreshold}(\beta, \gamma) := \begin{cases} \beta + \gamma &\text{if } \beta < -\gamma \\ 0 &\text{if } -\gamma \leq \beta \leq \gamma \\ \beta - \gamma &\text{if } \beta > \gamma. \end{cases}$$

Esta regla de actualización tiene las siguientes dos propiedades inspiradoras, que se explican a continuación:

1. En el caso de limitación $r_{\text{L1}} \to 0$ (es decir, sin regularización L1), esta regla de actualización es idéntica a la regla de actualización original.

2. Esta regla de actualización se puede interpretar como la aplicación de un operador de proximidad cuyo punto fijo es la solución al problema de minimización regularizado L1.

$$\underset{\beta - \beta^{(t)} \in \text{span}\{ \text{onehot}(j^{(t)}) \}}{\text{arg min}} \left( -\ell(\beta \,;\, \mathbf{x}, \mathbf{y}) + r_{\text{L1}} \left\lVert \beta \right\rVert_1 \right).$$

El caso degenerado $r_{\text{L1}} = 0$ recupera la regla de actualización original

Para ver (1), tenga en cuenta que si $r_{\text{L1}} = 0$ entonces $\gamma^{(t)} = 0$, por lo tanto

$$\begin{align*} \left(\beta_{\text{reg}}^{(t+1)}\right)_{j^{(t)}} &= \text{SoftThreshold} \left( \beta^{(t)}_{j^{(t)}} - \alpha\, u^{(t)} ,\ 0 \right) \\ &= \beta^{(t)}_{j^{(t)}} - \alpha\, u^{(t)}. \end{align*}$$

Por eso

$$\begin{align*} \beta_{\text{reg}}^{(t+1)} &= \beta^{(t)} - \alpha\, u^{(t)} \,\text{onehot}(j^{(t)}) \\ &= \beta^{(t+1)}. \end{align*}$$

Operador de proximidad cuyo punto fijo es el MLE regularizado

Para ver (2), primero tenga en cuenta (ver Wikipedia) que para cualquier $\gamma > 0$, la siguiente regla de actualización

$$\left(\beta_{\text{exact-prox}, \gamma}^{(t+1)}\right)_{j^{(t)}} := \text{prox}_{\gamma \lVert \cdot \rVert_1} \left( \beta^{(t)}_{j^{(t)}} + \frac{\gamma}{r_{\text{L1}}} \left( \left( \nabla_\beta\, \ell(\beta \,;\, \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)}} \right)_{j^{(t)}} \right)$$

satisface (2), donde $\text{prox}$ es el operador de proximidad (ver Yu, donde este operador se denota $\mathsf{P}$). El lado derecho de la ecuación anterior se calcula aquí:

$$\left(\beta_{\text{exact-prox}, \gamma}^{(t+1)}\right)_{j^{(t)}} = \text{SoftThreshold} \left( \beta^{(t)}_{j^{(t)}} + \frac{\gamma}{r_{\text{L1}}} \left( \left( \nabla_\beta\, \ell(\beta \,;\, \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)}} \right)_{j^{(t)}} ,\ \gamma \right).$$

En particular, al establecer $\gamma = \gamma^{(t)} = -\frac{\alpha\, r_{\text{L1}}}{\left(H^{(t)}\right)_{j^{(t)}, j^{(t)}}}$ (tenga en cuenta que $\gamma^{(t)} > 0$ siempre que la probabilidad logarítmica negativa sea convexa), obtenemos la regla de actualización

$$\left(\beta_{\text{exact-prox}, \gamma^{(t)}}^{(t+1)}\right)_{j^{(t)}} = \text{SoftThreshold} \left( \beta^{(t)}_{j^{(t)}} - \alpha \frac{ \left( \left( \nabla_\beta\, \ell(\beta \,;\, \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)}} \right)_{j^{(t)}} }{ \left(H^{(t)}\right)_{j^{(t)}, j^{(t)}} } ,\ \gamma^{(t)} \right).$$

Luego reemplazamos el gradiente exacto $\left( \nabla_\beta\, \ell(\beta \,;\, \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)}}$ con su aproximación $s^{(t)}$, y obtenemos

$$\begin{align*} \left(\beta_{\text{exact-prox}, \gamma^{(t)}}^{(t+1)}\right)_{j^{(t)}} &\approx \text{SoftThreshold} \left( \beta^{(t)}_{j^{(t)}} - \alpha \frac{ \left(s^{(t)}\right)_{j^{(t)}} }{ \left(H^{(t)}\right)_{j^{(t)}, j^{(t)}} } ,\ \gamma^{(t)} \right) \\ &= \text{SoftThreshold} \left( \beta^{(t)}_{j^{(t)}} - \alpha\, u^{(t)} ,\ \gamma^{(t)} \right). \end{align*}$$

Por eso

$$\beta_{\text{exact-prox}, \gamma^{(t)}}^{(t+1)} \approx \beta_{\text{reg}}^{(t+1)}.$$

Derivación de hechos del GLM

En esta sección, declaramos con todo detalle y derivamos los resultados sobre los GLM que se utilizan en las secciones anteriores. Luego, utilizamos gradients de TensorFlow para verificar numéricamente las fórmulas derivadas para el gradiente de la probabilidad logarítmica y la información de Fisher.

Puntuación e información de Fisher

Considere una familia de distribuciones de probabilidad parametrizadas por el vector de parámetros $\theta$, que tiene densidades de probabilidad $\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T}}$. La puntuación de un resultado $y$ en el vector de parámetros $\theta_0$ se define como el gradiente de la probabilidad logarítmica de $y$ (evaluada en $\theta_0$), es decir,

$$\text{score}(y, \theta_0) := \left[\nabla_\theta\, \log p(y | \theta)\right]_{\theta=\theta_0}.$$

Afirmación: La expectativa de la puntuación es cero

En condiciones de regularidad leves (que nos permiten pasar la diferenciación bajo la integral),

$$\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[\text{score}(Y, \theta_0)\right] = 0.$$

Prueba

Tenemos

$$\begin{align*} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[\text{score}(Y, \theta_0)\right] &:=\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[\left(\nabla_\theta \log p(Y|\theta)\right)_{\theta=\theta_0}\right] \\ &\stackrel{\text{(1)}}{=} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[\frac{\left(\nabla_\theta p(Y|\theta)\right)_{\theta=\theta_0}}{p(Y|\theta=\theta_0)}\right] \\ &\stackrel{\text{(2)}}{=} \int_{\mathcal{Y}} \left[\frac{\left(\nabla_\theta p(y|\theta)\right)_{\theta=\theta_0}}{p(y|\theta=\theta_0)}\right] p(y | \theta=\theta_0)\, dy \\ &= \int_{\mathcal{Y}} \left(\nabla_\theta p(y|\theta)\right)_{\theta=\theta_0}\, dy \\ &\stackrel{\text{(3)}}{=} \left[\nabla_\theta \left(\int_{\mathcal{Y}} p(y|\theta)\, dy\right) \right]_{\theta=\theta_0} \\ &\stackrel{\text{(4)}}{=} \left[\nabla_\theta\, 1 \right]_{\theta=\theta_0} \\ &= 0, \end{align*}$$

donde hemos utilizado: (1) regla de la cadena para la diferenciación, (2) definición de expectativa, (3) paso de la diferenciación bajo el signo integral (mediante el uso de las condiciones de regularidad), (4) la integral de una densidad de probabilidad es 1.

Afirmación (información de Fisher): La varianza de la puntuación es igual a la probabilidad hessiana negativa esperada de la probabilidad logarítmica.

En condiciones de regularidad leves (que nos permiten pasar la diferenciación bajo la integral),

$$\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \text{score}(Y, \theta_0) \text{score}(Y, \theta_0)^\top \right] = -\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \left(\nabla_\theta^2 \log p(Y | \theta)\right)_{\theta=\theta_0} \right]$$

donde $\nabla_\theta^2 F$ denota la matriz hessiana, cuya entrada $(i, j)$ es $\frac{\partial^2 F}{\partial \theta_i \partial \theta_j}$.

El lado izquierdo de esta ecuación se llama información de Fisher de la familia $\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T}}$ en el vector de parámetros $\theta_0$.

Prueba de la afirmación

Tenemos

$$\begin{align*} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \left(\nabla_\theta^2 \log p(Y | \theta)\right)_{\theta=\theta_0} \right] &\stackrel{\text{(1)}}{=} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \left(\nabla_\theta^\top \frac{ \nabla_\theta p(Y | \theta) }{ p(Y|\theta) }\right)_{\theta=\theta_0} \right] \\ &\stackrel{\text{(2)}}{=} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \frac{ \left(\nabla^2_\theta p(Y | \theta)\right)_{\theta=\theta_0} }{ p(Y|\theta=\theta_0) } - \left(\frac{ \left(\nabla_\theta\, p(Y|\theta)\right)_{\theta=\theta_0} }{ p(Y|\theta=\theta_0) }\right) \left(\frac{ \left(\nabla_\theta\, p(Y|\theta)\right)_{\theta=\theta_0} }{ p(Y|\theta=\theta_0) }\right)^\top \right] \\ &\stackrel{\text{(3)}}{=} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \frac{ \left(\nabla^2_\theta p(Y | \theta)\right)_{\theta=\theta_0} }{ p(Y|\theta=\theta_0) } - \text{score}(Y, \theta_0) \,\text{score}(Y, \theta_0)^\top \right], \end{align*}$$

donde hemos usado (1) la regla de la cadena para la diferenciación, (2) la regla del cociente para la diferenciación, (3) la regla de la cadena nuevamente, a la inversa.

Para completar la demostración, basta demostrar lo siguiente:

$$\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \frac{ \left(\nabla^2_\theta p(Y | \theta)\right)_{\theta=\theta_0} }{ p(Y|\theta=\theta_0) } \right] \stackrel{\text{?}}{=} 0.$$

Para hacer eso, pasamos la diferenciación bajo el signo integral dos veces:

$$\begin{align*} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \frac{ \left(\nabla^2_\theta p(Y | \theta)\right)_{\theta=\theta_0} }{ p(Y|\theta=\theta_0) } \right] &= \int_{\mathcal{Y}} \left[ \frac{ \left(\nabla^2_\theta p(y | \theta)\right)_{\theta=\theta_0} }{ p(y|\theta=\theta_0) } \right] \, p(y | \theta=\theta_0)\, dy \\ &= \int_{\mathcal{Y}} \left(\nabla^2_\theta p(y | \theta)\right)_{\theta=\theta_0} \, dy \\ &= \left[ \nabla_\theta^2 \left( \int_{\mathcal{Y}} p(y | \theta) \, dy \right) \right]_{\theta=\theta_0} \\ &= \left[ \nabla_\theta^2 \, 1 \right]_{\theta=\theta_0} \\ &= 0. \end{align*}$$

Lema sobre la derivada de la función de partición logarítmica

Si $a$, $b$ y $c$ son funciones con valores escalares, $c$ dos veces diferenciables, de modo que la familia de distribuciones $\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T}}$ definido por

$$p(y|\theta) = a(y) \exp\left(b(y)\, \theta - c(\theta)\right)$$

cumple las condiciones de regularidad leves que permiten pasar la diferenciación con respecto a $\theta$ bajo una integral con respecto a $y$, entonces

$$\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ b(Y) \right] = c'(\theta_0)$$

y

$$\text{Var}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ b(Y) \right] = c''(\theta_0).$$

(Aquí $'$ denota diferenciación, por lo que $c'$ y $c''$ son la primera y la segunda derivada de $c$).

Prueba

Para esta familia de distribuciones, tenemos $\text{score}(y, \theta_0) = b(y) - c'(\theta_0)$. La primera ecuación se deriva entonces del hecho $\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ \text{score}(y, \theta_0) \right] = 0$. A continuación, tenemos

$$\begin{align*} \text{Var}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ b(Y) \right] &= \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ \left(b(Y) - c'(\theta_0)\right)^2 \right] \\ &= \text{the one entry of } \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ \text{score}(y, \theta_0) \text{score}(y, \theta_0)^\top \right] \\ &= \text{the one entry of } -\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ \left(\nabla_\theta^2 \log p(\cdot | \theta)\right)_{\theta=\theta_0} \right] \\ &= -\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ -c''(\theta_0) \right] \\ &= c''(\theta_0). \end{align*}$$

Familia exponencial de distribuciones sobredispersas

Una familia exponencial (escalar) de distribuciones sobredispersas es una familia de distribuciones cuyas densidades toman la forma

$$p_{\text{OEF}(m, T)}(y\, |\, \theta, \phi) = m(y, \phi) \exp\left(\frac{\theta\, T(y) - A(\theta)}{\phi}\right),$$

donde $m$ y $T$ son funciones escalares conocidas, y $\theta$ y $\phi$ son parámetros escalares.

[Tenga en cuenta que $A$ está sobredeterminado: para cualquier $\phi_0$, la función $A$ está completamente determinada por la restricción de que $\int p_{\text{OEF}(m, T)}(y\ |\ \theta, \phi=\phi_0)\, dy = 1$ para todos los $\theta$. Los $A$ producidos por diferentes valores de $\phi_0$ deben ser todos iguales, lo que impone una restricción a las funciones $m$ y $T$].

Media y varianza del estadístico suficiente

En las mismas condiciones que en el "Lema sobre la derivada de la función de partición logarítmica", tenemos

$$\mathbb{E}_{Y \sim p_{\text{OEF}(m, T)}(\cdot | \theta, \phi)} \left[ T(Y) \right] = A'(\theta)$$

y

$$\text{Var}_{Y \sim p_{\text{OEF}(m, T)}(\cdot | \theta, \phi)} \left[ T(Y) \right] = \phi A''(\theta).$$

Prueba

De acuerdo con el "Lema sobre la derivada de la función de partición logarítmica", tenemos

$$\mathbb{E}_{Y \sim p_{\text{OEF}(m, T)}(\cdot | \theta, \phi)} \left[ \frac{T(Y)}{\phi} \right] = \frac{A'(\theta)}{\phi}$$

y

$$\text{Var}_{Y \sim p_{\text{OEF}(m, T)}(\cdot | \theta, \phi)} \left[ \frac{T(Y)}{\phi} \right] = \frac{A''(\theta)}{\phi}.$$

El resultado se deriva entonces del hecho de que la expectativa es lineal ($\mathbb{E}[aX] = a\mathbb{E}[X]$) y la varianza es homogénea de grado 2 ($\text{Var}[aX] = a^2 \,\text{Var}[X]$).

Modelo lineal generalizado

En un modelo lineal generalizado, una distribución predictiva para la variable de respuesta $Y$ se asocia con un vector de predictores observados $x$. La distribución es miembro de una familia exponencial de distribuciones sobredispersas y el parámetro $\theta$ se reemplaza por $h(\eta)$ donde $h$ es una función conocida, $\eta := x^\top \beta$ es la llamada respuesta lineal y $\beta$ es un vector de parámetros (coeficientes de regresión) que se deben aprender. En general, el parámetro de dispersión $\phi$ también se puede aprender, pero en nuestra configuración trataremos $\phi$ como si fuera conocido. Entonces nuestra configuración será la siguiente:

$$Y \sim p_{\text{OEF}(m, T)}(\cdot\, |\, \theta = h(\eta), \phi)$$

donde la estructura del modelo se caracteriza por la distribución $p_{\text{OEF}(m, T)}$ y la función $h$ que convierte la respuesta lineal en parámetros.

Tradicionalmente, la asignación de la respuesta lineal $\eta$ al significado $\mu := \mathbb{E}_{Y \sim p_{\text{OEF}(m, T)}(\cdot\, |\, \theta = h(\eta), \phi)}\left[ Y\right]$ se denota así:

$$\mu = g^{-1}(\eta).$$

Se requiere que esta asignación sea uno a uno, y su inverso, $g$, se denomina función de enlace para este GLM. Normalmente, se describe un GLM nombrando su función de enlace y su familia de distribuciones; por ejemplo, un "GLM con distribución de Bernoulli y función de enlace logit" (también conocido como modelo de regresión logística). Para caracterizar completamente el GLM, también se debe especificar la función $h$. Si $h$ es la identidad, entonces se dice que $g$ es la función de enlace canónico.

Afirmación: $h'$ se expresa en términos del estadístico suficiente

Defina

$${\text{Mean}_T}(\eta) := \mathbb{E}_{Y \sim p_{\text{OEF}(m, T)}(\cdot | \theta = h(\eta), \phi)} \left[ T(Y) \right]$$

y

$${\text{Var}_T}(\eta) := \text{Var}_{Y \sim p_{\text{OEF}(m, T)}(\cdot | \theta = h(\eta), \phi)} \left[ T(Y) \right].$$

Entonces, tenemos

$$h'(\eta) = \frac{\phi\, {\text{Mean}_T}'(\eta)}{{\text{Var}_T}(\eta)}.$$

Prueba

Mediante la "Media y varianza del estadístico suficiente", tenemos

$${\text{Mean}_T}(\eta) = A'(h(\eta)).$$

Al derivar con la regla de la cadena obtenemos $${\text{Mean}_T}'(\eta) = A''(h(\eta))\, h'(\eta),$$

y mediante la "Media y varianza del estadístico suficiente",

$$\cdots = \frac{1}{\phi} {\text{Var}_T}(\eta)\ h'(\eta).$$

La conclusión es la siguiente.

Ajuste de parámetros del GLM a los datos

Las propiedades derivadas anteriormente se prestan muy bien para ajustar los parámetros del GLM $\beta$ a un conjunto de datos. Los métodos cuasi-Newton, como la puntuación de Fisher, se basan en el gradiente del logaritmo de probabilidad y la información de Fisher, que ahora mostramos se puede calcular de manera especialmente eficiente para un GLM.

Supongamos que hemos observado vectores predictores $x_i$ y respuestas escalares asociadas $y_i$. En forma matricial, diremos que hemos observado predictores $\mathbf{x}$ y respuesta $\mathbf{y}$, donde $\mathbf{x}$ es la matriz cuya $i$ésima fila es $x_i^ \top$ y $\mathbf{y}$ es el vector cuyo $i$ésimo elemento es $y_i$. La probabilidad logarítmica de los parámetros $\beta$ es entonces la siguiente:

$$\ell(\beta\, ;\, \mathbf{x}, \mathbf{y}) = \sum_{i=1}^{N} \log p_{\text{OEF}(m, T)}(y_i\, |\, \theta = h(x_i^\top \beta), \phi).$$

Para una sola muestra de datos

Para simplificar la notación, consideremos primero el caso de un único punto de datos, $N=1$; luego, extenderemos al caso general por aditividad.

Gradiente

Tenemos

$$\begin{align*} \ell(\beta\, ;\, x, y) &= \log p_{\text{OEF}(m, T)}(y\, |\, \theta = h(x^\top \beta), \phi) \\ &= \log m(y, \phi) + \frac{\theta\, T(y) - A(\theta)}{\phi}, \quad\text{where}\ \theta = h(x^\top \beta). \end{align*}$$

Por lo tanto, según la regla de la cadena,

$$\nabla_\beta \ell(\beta\, ; \, x, y) = \frac{T(y) - A'(\theta)}{\phi}\, h'(x^\top \beta)\, x.$$

Por separado, mediante la "Media y varianza del estadístico suficiente", tenemos $A'(\theta) = {\text{Media}_T}(x^\top \beta)$. Por lo tanto, según "Afirmación: $h'$ ser expresa en términos del estadístico suficiente", tenemos

$$\cdots = \left(T(y) - {\text{Mean}_T}(x^\top \beta)\right) \frac{{\text{Mean}_T}'(x^\top \beta)}{{\text{Var}_T}(x^\top \beta)} \,x.$$

Hessiano

Al derivar por segunda vez, de acuerdo con la regla del producto obtenemos

$$\begin{align*} \nabla_\beta^2 \ell(\beta\, ;\, x, y) &= \left[ -A''(h(x^\top \beta))\, h'(x^\top \beta) \right] h'(x^\top \beta)\, x x^\top + \left[ T(y) - A'(h(x^\top \beta)) \right] h''(x^\top \beta)\, xx^\top ] \\ &= \left( -{\text{Mean}_T}'(x^\top \beta)\, h'(x^\top \beta) + \left[T(y) - A'(h(x^\top \beta))\right] \right)\, x x^\top. \end{align*}$$

Información de Fisher

Mediante la "Media y varianza del estadístico suficiente", tenemos

$$\mathbb{E}_{Y \sim p_{\text{OEF}(m, T)}(\cdot | \theta = h(x^\top \beta), \phi)} \left[ T(y) - A'(h(x^\top \beta)) \right] = 0.$$

Por eso

$$\begin{align*} \mathbb{E}_{Y \sim p_{\text{OEF}(m, T)}(\cdot | \theta = h(x^\top \beta), \phi)} \left[ \nabla_\beta^2 \ell(\beta\, ;\, x, y) \right] &= -{\text{Mean}_T}'(x^\top \beta)\, h'(x^\top \beta) x x^\top \\ &= -\frac{\phi\, {\text{Mean}_T}'(x^\top \beta)^2}{{\text{Var}_T}(x^\top \beta)}\, x x^\top. \end{align*}$$

Para múltiples muestras de datos

Ahora extendemos el caso $N=1$ al caso general. Digamos que $\boldsymbol{\eta} := \mathbf{x} \beta$ denota el vector cuya $i$ésima coordenada es la respuesta lineal de la $i$ésima muestra de datos. Supongamos que $\mathbf{T}$ (resp. ${\textbf{Mean}_T}$, resp. ${\textbf{Var}_T}$) denota la función transmitida (vectorizada) que aplica la función de valor escalar $T$ (resp. ${\text{Mean}_T}$, resp. ${\text{Var}_T}$) a cada coordenada. Entonces tenemos

$$\begin{align*} \nabla_\beta \ell(\beta\, ;\, \mathbf{x}, \mathbf{y}) &= \sum_{i=1}^{N} \nabla_\beta \ell(\beta\, ;\, x_i, y_i) \\ &= \sum_{i=1}^{N} \left(T(y) - {\text{Mean}_T}(x_i^\top \beta)\right) \frac{{\text{Mean}_T}'(x_i^\top \beta)}{{\text{Var}_T}(x_i^\top \beta)} \, x_i \\ &= \mathbf{x}^\top \,\text{diag}\left(\frac{ {\textbf{Mean}_T}'(\mathbf{x} \beta) }{ {\textbf{Var}_T}(\mathbf{x} \beta) }\right) \left(\mathbf{T}(\mathbf{y}) - {\textbf{Mean}_T}(\mathbf{x} \beta)\right) \\ \end{align*}$$

y

$$\begin{align*} \mathbb{E}_{Y_i \sim p_{\text{OEF}(m, T)}(\cdot | \theta = h(x_i^\top \beta), \phi)} \left[ \nabla_\beta^2 \ell(\beta\, ;\, \mathbf{x}, \mathbf{Y}) \right] &= \sum_{i=1}^{N} \mathbb{E}_{Y_i \sim p_{\text{OEF}(m, T)}(\cdot | \theta = h(x_i^\top \beta), \phi)} \left[ \nabla_\beta^2 \ell(\beta\, ;\, x_i, Y_i) \right] \\ &= \sum_{i=1}^{N} -\frac{\phi\, {\text{Mean}_T}'(x_i^\top \beta)^2}{{\text{Var}_T}(x_i^\top \beta)}\, x_i x_i^\top \\ &= -\mathbf{x}^\top \,\text{diag}\left( \frac{ \phi\, {\textbf{Mean}_T}'(\mathbf{x} \beta)^2 }{ {\textbf{Var}_T}(\mathbf{x} \beta) }\right)\, \mathbf{x}, \end{align*}$$

donde las fracciones denotan división por elementos.

Cómo verificar las fórmulas numéricamente

Ahora usamos tf.gradients para verificar numéricamente la fórmula anterior para el gradiente de la probabilidad logarítmica y verificamos la fórmula para la información de Fisher con una estimación de Monte Carlo con ayuda de tf.hessians:

def VerifyGradientAndFIM():
  model = tfp.glm.BernoulliNormalCDF()
  model_matrix = np.array([[1., 5, -2],
                           [8, -1, 8]])

  def _naive_grad_and_hessian_loss_fn(x, response):
    # Computes gradient and Hessian of negative log likelihood using autodiff.
    predicted_linear_response = tf.linalg.matvec(model_matrix, x)
    log_probs = model.log_prob(response, predicted_linear_response)
    grad_loss = tf.gradients(-log_probs, [x])[0]
    hessian_loss = tf.hessians(-log_probs, [x])[0]
    return [grad_loss, hessian_loss]

  def _grad_neg_log_likelihood_and_fim_fn(x, response):
    # Computes gradient of negative log likelihood and Fisher information matrix
    # using the formulas above.
    predicted_linear_response = tf.linalg.matvec(model_matrix, x)
    mean, variance, grad_mean = model(predicted_linear_response)

    v = (response - mean) * grad_mean / variance
    grad_log_likelihood = tf.linalg.matvec(model_matrix, v, adjoint_a=True)
    w = grad_mean**2 / variance

    fisher_info = tf.linalg.matmul(
        model_matrix,
        w[..., tf.newaxis] * model_matrix,
        adjoint_a=True)
    return [-grad_log_likelihood, fisher_info]

  @tf.function(autograph=False)
  def compute_grad_hessian_estimates():
    # Monte Carlo estimate of E[Hessian(-LogLikelihood)], where the expectation is
    # as written in "Claim (Fisher information)" above.
    num_trials = 20
    trial_outputs = []
    np.random.seed(10)
    model_coefficients_ = np.random.random(size=(model_matrix.shape[1],))
    model_coefficients = tf.convert_to_tensor(model_coefficients_)
    for _ in range(num_trials):
      # Sample from the distribution of `model`
      response = np.random.binomial(
          1,
          scipy.stats.norm().cdf(np.matmul(model_matrix, model_coefficients_))
      ).astype(np.float64)
      trial_outputs.append(
          list(_naive_grad_and_hessian_loss_fn(model_coefficients, response)) +
          list(
              _grad_neg_log_likelihood_and_fim_fn(model_coefficients, response))
      )

    naive_grads = tf.stack(
        list(naive_grad for [naive_grad, _, _, _] in trial_outputs), axis=0)
    fancy_grads = tf.stack(
        list(fancy_grad for [_, _, fancy_grad, _] in trial_outputs), axis=0)

    average_hess = tf.reduce_mean(tf.stack(
        list(hess for [_, hess, _, _] in trial_outputs), axis=0), axis=0)
    [_, _, _, fisher_info] = trial_outputs[0]
    return naive_grads, fancy_grads, average_hess, fisher_info
  
  naive_grads, fancy_grads, average_hess, fisher_info = [
      t.numpy() for t in compute_grad_hessian_estimates()]

  print("Coordinatewise relative error between naively computed gradients and"
        " formula-based gradients (should be zero):\n{}\n".format(
            (naive_grads - fancy_grads) / naive_grads))

  print("Coordinatewise relative error between average of naively computed"
        " Hessian and formula-based FIM (should approach zero as num_trials"
        " -> infinity):\n{}\n".format(
                (average_hess - fisher_info) / average_hess))
    
VerifyGradientAndFIM()

Coordinatewise relative error between naively computed gradients and formula-based gradients (should be zero):
[[2.08845965e-16 1.67076772e-16 2.08845965e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [2.08845965e-16 1.67076772e-16 2.08845965e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [2.08845965e-16 1.67076772e-16 2.08845965e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [2.08845965e-16 1.67076772e-16 2.08845965e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [2.08845965e-16 1.67076772e-16 2.08845965e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [2.08845965e-16 1.67076772e-16 2.08845965e-16]
 [1.96118673e-16 3.13789877e-16 1.96118673e-16]
 [2.08845965e-16 1.67076772e-16 2.08845965e-16]]

Coordinatewise relative error between average of naively computed Hessian and formula-based FIM (should approach zero as num_trials -> infinity):
[[0.00072369 0.00072369 0.00072369]
 [0.00072369 0.00072369 0.00072369]
 [0.00072369 0.00072369 0.00072369]]

Referencias

[1]: Guo-Xun Yuan, Chia-Hua Ho y Chih-Jen Lin. An Improved GLMNET for L1-regularized Logistic Regression. Journal of Machine Learning Research, 13, 2012. http://www.jmlr.org/papers/volume13/yuan12a/yuan12a.pdf

[2]: skd. Derivation of Soft Thresholding Operator. 2018. https://math.stackexchange.com/q/511106

[3]: Contribuyentes de Wikipedia. Proximal gradient methods for learning. Wikipedia, The Free Encyclopedia, 2018. https://en.wikipedia.org/wiki/Proximal_gradient_methods_for_learning

[4]: Yao-Liang Yu. The Proximity Operator. https://www.cs.cmu.edu/~suvrit/teach/yaoliang_proximity.pdf