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 lineares generalizados

Veja em TensorFlow.org

Executar no Google Colab

Ver fonte em GitHub

Baixar notebook

Neste notebook, apresentamos os modelos lineares generalizados usando um exemplo funcional, que será resolvido de duas formas diferentes usando dois algoritmos para ajustar de forma eficiente modelos lineares generalizados no TensorFlow Probability: pontuação de Fisher para dados densos e método do gradiente proximal com coordenadas. Comparamos os coeficientes ajustados aos coeficientes verdadeiros e, no caso do método do gradiente proximal com coordenadas, à saída do algoritmo glmnet similar do R. Por fim, fornecemos maiores detalhes matemáticos e derivações de diversas propriedades essenciais dos modelos lineares generalizados.

Contexto

Um modelo linear generalizado (GLM, na sigla em inglês) é um modelo linear ($\eta = x^\top \beta$) encapsulado em uma transformação (função de ligação) que conta com uma distribuição de respostas a partir de uma família de exponenciais. A escolha da função de ligação e da distribuição de respostas é bastante flexível, proporcionando grande capacidade de expressão para os modelos lineares generalizados. Confira os detalhes completos, incluindo uma apresentação sequencial de todas as definições e resultados até a construção dos modelos lineares generalizados em notação não ambígua, na seção "Derivação de fatos sobre modelos lineares generalizados" abaixo. Vamos resumir:

Em um modelo linear generalizado, uma distribuição preditiva para a variável de resposta $Y$ está associada a um vetor de preditores observados $x$. A distribuição tem a seguinte 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*}$$

Aqui, $\beta$ são os parâmetros ("pesos"), $\phi$ é um hiperparâmetro que representa a dispersão ("variância"), e $m$, $h$, $T$, $A$ são caracterizados pela família de modelos especificados pelo usuário.

A média de $Y$ depende de $x$ pela composição da resposta linear $\eta$ e pela função de ligação (inversa), isto é:

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

em que $g$ é a função de ligação. No TFP, a função de ligação e a família de modelos são especificadas conjuntamente por uma subclasse de tfp.glm.ExponentialFamily. Confira alguns exemplos:

• tfp.glm.Normal, também chamada de "regressão linear"

• tfp.glm.Bernoulli, também chamada de "regressão logística"

• tfp.glm.Poisson, também chamada de "regressão de Poisson"

• tfp.glm.BernoulliNormalCDF, também chamada de "regressão probit"

O TFP prefere nomear as famílias de modelos de acordo com a distribuição ao longo de Y em vez da função de ligação, já que as distribuições tfp.Distribution são funções de primeira classe. Se o nome da subclasse de tfp.glm.ExponentialFamily contiver uma segunda palavra, isso indica uma função de ligação não canônica.

Os modelos lineares generalizados têm propriedades marcantes que permitem uma implementação eficiente do estimador de máxima verossimilhança. Uma dessas propriedades são fórmulas simples para o gradiente da log-verossimilhança $\ell$ e para a matriz de informação de Fisher, que é o valor esperado da Hessiana da log-verossimilhança negativa ao fazer uma reamostragem da resposta para os mesmos preditores, isto é:

$$\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*}$$

em que $\mathbf{x}$ é a matriz cuja $i$ésima linha é o vetor de preditores para a $i$ésima amostra de dados, e $\mathbf{y}$ é um vetor cuja $i$ésima coordenada é a resposta observada para a $i$ésima amostra de dados. Aqui (grosso modo), ${\text{Mean}_T}(\eta) := \mathbb{E}[T(Y),|,\eta]$ e ${\text{Var}_T}(\eta) := \text{Var}[T(Y),|,\eta]$, e o negrito denota a vetorização dessas funções. Confira os detalhes completos sobre as expectativas e variâncias dessas distribuições na seção "Derivação de fatos sobre modelos lineares generalizados" abaixo.

Exemplo

Nesta seção, vamos descrever e demonstrar brevemente dois algoritmos integrados de ajuste de modelos lineares generalizados no TensorFlow Probability: pontuação de Fisher (tfp.glm.fit) e método do gradiente proximal com coordenadas (tfp.glm.fit_sparse).

Dataset sintético

Vamos fingir o carregamento de um dataset de treinamento.

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

Observação: estabeleça conexão a um runtime local.

Neste notebook, compartilhamos dados entre os kernels do Python e do R usando arquivos locais. Para permitir esse compartilhamento, use runtimes na mesma máquina na qual você tenha permissão de ler e gravar arquivos locais.

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=',')

Sem regularização L1

A função tfp.glm.fit implementa a pontuação de Fisher, que recebe como alguns de seus argumentos:

• model_matrix = $\mathbf{x}$

• response = $\mathbf{y}$

• model = callable que, dado o argumento $\boldsymbol{\eta}$, retorna a tupla $\left( {\textbf{Mean}_T}(\boldsymbol{\eta}), {\textbf{Var}_T}(\boldsymbol{\eta}), {\textbf{Mean}_T}'(\boldsymbol{\eta}) \right)$.

Recomendamos que model seja uma instância da classe tfp.glm.ExponentialFamily. Existem diversas implementações pré-criadas disponíveis, então, para a maioria dos modelos lineares generalizados comuns, não é necessário escrever 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

Detalhes matemáticos

A pontuação de Fisher é uma modificação do método de Newton para encontrar a estimativa da máxima verossimilhança:

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

O método padrão de Newton, que é a procura por zeros do gradiente da log-verossimilhança, seguiria a regra de atualização:

$$\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)}}$$

em que $\alpha \in (0, 1]$ é uma taxa de aprendizado usada para controlar o tamanho do passo.

Na pontuação de Fisher, substituímos a Hessiana pela matriz de informação de Fisher negativa:

$$\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*}$$

[Observe que, aqui, $\mathbf{Y} = (Y_i)_{i=1}^{n}$ é aleatório, enquanto $\mathbf{y}$ ainda é o vetor de respostas observadas.]

Segundo as fórmulas na seção "Ajuste dos parâmetros de modelos lineares generalizados de acordo com os dados" abaixo, isso é simplificado para:

$$\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*}$$

Com regularização L1

tfp.glm.fit_sparse implementa um ajustador de modelos lineares generalizados mais adequado para datasets esparsos baseado no algoritmo em Yuan, Ho e Lin, 2012. Confira alguns recursos:

• Regularização L1

• Sem inversão de matrizes

• Poucas avaliações do gradiente e da Hessiana

Primeiro, vamos apresentar um exemplo de uso do código. Os detalhes do algoritmo são mais bem explicados na seção "Detalhes do algoritmo para tfp.glm.fit_sparse" abaixo.

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:

Observe que os coeficientes aprendidos têm o mesmo padrão de esparsidade que os coeficientes verdadeiros.

# 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=',')

Comparação com glmnet do R

Vamos comparar a saída do método do gradiente proximal com coordenadas ao glmnet do R, que usa um algoritmo similar.

OBSERVAÇÃO: para executar o código nesta seção, você precisa mudar para um runtime de Colab do 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)

Comparação entre R, TFP e coeficientes verdadeiros (observação: volte para o kernel do 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

Detalhes do algoritmo para tfp.glm.fit_sparse

Apresentamos o algoritmo como uma sequência de três modificações do método de Newton. Em cada uma delas, a regra de atualização de $\beta$ é baseada em um vetor $s$ e em uma matriz $H$, que faz a aproximação do gradiente e da Hessiana da log-verossimilhança. No passo $t$, escolhemos uma coordenada $j^{(t)}$ a ser alterada e atualizamos $\beta$ de acordo com a regra de atualização:

$$\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*}$$

Essa atualização é um passo similar ao Newton com taxa de aprendizado $\alpha$. Exceto pela parte final (regularização L1), as modificações abaixo diferem apenas na forma como $s$ e $H$ são atualizados.

Ponto de partida: método de Newton com coordenadas

No método de Newton com coordenadas, definimos $s$ e $H$ como o gradiente verdadeiro e a Hessiana da log-verossimilhança:

$$\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 avaliações do gradiente e da Hessiana

Geralmente, é caro computar o gradiente e a Hessiana da log-verossimilhança, então costuma valer a pena fazer a aproximação deles da seguinte forma:

• Geralmente, fazemos a aproximação da Hessiana como constante localmente e fazemos a aproximação do gradiente para a primeira ordem usando a Hessiana (aproximada):

$$\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*}$$

• Ocasionalmente, faça um passo de atualização padrão conforme acima, definindo $s^{(t+1)}$ como o gradiente exato e $H^{(t+1)}$ como a Hessiana exata da log-verossimilhança, avaliada em $\beta^{(t+1)}$.

Substituta a informação de Fisher negativa pela Hessiana

Para reduzir ainda mais o custo dos passos de atualização padrão, podemos definir $H$ como a matriz de informação de Fisher negativa (eficiente do ponto de vista computacional ao usar as fórmulas da seção "Ajuste dos parâmetros de modelos lineares generalizados de acordo com os dados" abaixo) em vez da Hessiana exata:

$$\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*}$$

Regularização L1 via método do gradiente proximal

Para incorporar a regularização L1, substituímos a regra de atualização

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

pela regra de atualização mais geral

$$\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*}$$

em que $r_{ ext{L1}} > 0$ é uma constante fornecida (o coeficiente de regularização L1), e $\text{SoftThreshold}$ é o operador de limiarização suave (soft threshold), 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}$$

Essa regra de atualização tem as duas propriedades inspiradoras abaixo:

1. No caso limitante $r_{\text{L1}} \to 0$ (ou seja, sem regularização L1), essa regra de atualização é idêntica à original.

2. Essa regra de atualização pode ser interpretada como a aplicação de um operador de proximidade cujo ponto fixo é a solução do problema de minimização com regularização 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).$$

O caso degenerado $r_{\text{L1}} = 0$ recupera a regra de atualização original

Para entender (1), observe que, se $r_{\text{L1}} = 0$, então $\gamma^{(t)} = 0$, portanto:

$$\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*}$$

Portanto:

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

Operador de proximidade cujo ponto fixo é a estimativa de máxima verossimilhança regularizada

Para entender (2), primeiro observe (confira a Wikipedia) que, para qualquer $\gamma > 0$, a regra de atualização

$$\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)$$

satisfaz (2), em que $\text{prox}$ é o operador de proximidade (confira Yu, em que esse operador é denotado como $\mathsf{P}$). O lado direito da equação acima é computado aqui:

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

Especificamente, ao definir $\gamma = \gamma^{(t)} = -\frac{\alpha\, r_{\text{L1}}}{\left(H^{(t)}\right)_{j^{(t)}, j^{(t)}}}$ (observe que (note that $\gamma^{(t)} > 0$, desde que a log-verossimilhança negativa seja convexa), obtemos a regra de atualização:

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

Em seguida, substitua o gradiente exato $\left( \nabla_\beta\, \ell(\beta \,;\, \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)}}$ por sua aproximação $s^{(t)}$, obtendo:

$$\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*}$$

Portanto:

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

Derivação de fatos sobre modelos lineares generalizados

Nesta seção, apresentamos maiores detalhes e derivamos os resultados sobre modelos lineares generalizados que foram usados nas seções anteriores. Em seguida, usamos gradients do TensorFlow para verificar numericamente as fórmulas derivadas para o gradiente da log-verossimilhança e da informação de Fisher.

Pontuação e informação de Fisher

Considere uma família de distribuições de probabilidade parametrizadas pelo vetor de parâmetros $\theta$, com densidades de probabilidade $\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T}}$. A pontuação de um resultado $y$ no vetor de parâmetros $\theta_0$ é definida como o gradiente da log-verossimilhança de $y$ (avaliada em $\theta_0$), isto é:

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

Afirmação: a expectativa da pontuação é zero

Sob condições de regularidade fracas (permitindo passar a diferenciação sob a integral),

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

Prova

Temos:

$$\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*}$$

em que usamos: (1) regra da cadeia para diferenciação, (2) definição da expectativa, (3) passagem da diferenciação sob o sinal da integral (usando as condições de regularidade), (4) a integral de uma densidade de probabilidade é igual a 1.

Afirmação (informação de Fisher): a variância da pontuação é igual à Hessiana esperada negativa da log-verossimilhança

Sob condições de regularidade fracas (permitindo passar a diferenciação sob a 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]$$

em que $\nabla_\theta^2 F$ denota a matriz Hessiana, cuja entrada $(i, j)$ é $\frac{\partial^2 F}{\partial \theta_i \partial \theta_j}$.

O lado esquerdo dessa equação é chamado de informação de Fisher da família $\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T}}$ no vetor de parâmetros $\theta_0$.

Prova da afirmação

Temos:

$$\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*}$$

em que usamos: (1) regra da cadeia para diferenciação, (2) regra do quociente para diferenciação, (3) regra da cadeia novamente, mas reversa.

Para completar a prova, é suficiente demonstrar que:

$$\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 fazer isso, passamos a diferenciação sob o sinal de integral duas vezes:

$$\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 a derivada da função de partição logarítmica

Se $a$, $b$ e $c$ são funções de valor escalar, $c$ sendo diferenciável duas vezes, de tal forma que a família de distribuições $\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T}}$ definida por

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

satisfaça as condições de regularidade fracas que permitem passar a diferenciação com relação a $\theta$ sob uma integral com relação a $y$, então:

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

e

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

(Aqui, $'$ denota a diferenciação, então $c'$ e $c''$ são a primeira e segunda derivadas de $c$. )

Prova

Para essa família de distribuições, temos $\text{score}(y, \theta_0) = b(y) - c'(\theta_0)$. A primeira equação vem então do fato de que $\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ \text{score}(y, \theta_0) \right] = 0$. Em seguida, temos:

$$\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*}$$

Família de exponenciais superdispersadas

Uma família de exponenciais superdispersadas (escalares) é uma família cujas densidades assumem a forma de:

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

em que $m$ e $T$ são funções de valor escalar conhecidas, e $\theta$ e $\phi$ são parâmetros escalares.

[Observe que $A$ é sobredeterminada: para qualquer $\phi_0$, a função $A$ é completamente determinada pela restrição de que $\int p_{\text{OEF}(m, T)}(y\ |\ \theta, \phi=\phi_0)\, dy = 1$ para todos os $\theta$. Todos os $A$ gerados por diferentes valores de $\phi_0$ devem ser iguais, o que coloca uma restrição nas funções $m$ e $T$.]

Média e variância da estatística suficiente

Sob as mesmas condições que em "Lema sobre a derivada da função de partição logarítmica" acima, temos:

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

e

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

Prova

De acordo com "Lema sobre a derivada da função de partição logarítmica", temos:

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

e

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

O resultado vem então do fato de que a expectativa é linear ($\mathbb{E}[aX] = a\mathbb{E}[X]$), e a variância é homogênea de grau 2 ($\text{Var}[aX] = a^2 \,\text{Var}[X]$).

Modelo linear generalizado

Em um modelo linear generalizado, uma distribuição preditiva para a variável de resposta $Y$ está associada a um vetor de preditores observados $x$. A distribuição é membro de uma família de exponenciais superdispersas, e o parâmetro $\theta$ é substituído por $h(\eta)$, em que $h$ é uma função conhecida, $\eta := x^\top \beta$ é a resposta linear, e $\beta$ é um vetor de parâmetros (coeficientes de regressão) a serem aprendidos. Em geral, o parâmetro de dispersão $\phi$ também pode ser aprendido, mas, em nossa configuração, tratamos $\phi$ como conhecido. Então, temos:

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

em que a estrutura do modelo é caracterizada pela distribuição $p_{\text{OEF}(m, T)}$ e pela função $h$, que converte a resposta linear em parâmetros.

Tradicionalmente, o mapeamento da resposta linear $\eta$ para a média $\mu := \mathbb{E}_{Y \sim p_{\text{OEF}(m, T)}(\cdot\, |\, \theta = h(\eta), \phi)}\left[ Y\right]$ é denotado assim:

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

Esse mapeamento precisa ser um-para-um obrigatoriamente, e seu mapeamento inverso, $g$, é chamado de função de ligação para esse modelo linear generalizado. Tipicamente, descrevemos um modelo linear generalizado nomeando sua função de ligação e sua família de distribuições – por exemplo: um "modelo linear generalizado com distribuição de Bernoulli e função de ligação logit" (também conhecido como modelo de regressão logística). Para caracterizar totalmente o modelo linear generalizado, a função $h$ também precisa ser especificada. Se $h$ é a identidade, então dizemos que $g$ é a função de ligação canônica.

Afirmação: expressando $h'$ em termos da estatística suficiente

Definir

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

e

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

Então, temos:

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

Prova

De acordo com "Média e variância da estatística suficiente", temos:

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

Fazendo a diferenciação com a regra da cadeia, obtemos $${\text{Mean}_T}'(\eta) = A''(h(\eta))\, h'(\eta),$$

e, de acordo com "Média e variância da estatística suficiente", temos:

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

A conclusão vem daí.

Ajuste dos parâmetros de modelos lineares generalizados de acordo com os dados

As propriedades derivadas acima são excelentes para ajustar os parâmetros $\beta$ de modelos lineares generalizados de acordo com um dataset. Os métodos quase-Newton, como a pontuação de Fisher, dependem do gradiente da log-verossimilhança e da informação de Fisher, que agora podemos computar de uma forma particularmente eficiente para um modelo linear generalizado.

Digamos que tenhamos os vetores de preditores observados $x_i$ e as respostas escalares associadas $y_i$. Na forma de matriz, vamos dizer que observamos os preditores $\mathbf{x}$ e a resposta $\mathbf{y}$, em que $\mathbf{x}$ é a matriz cuja $i$ésima linha é $x_i^\top$, e $\mathbf{y}$ é o vetor cujo $i$ésimo elemento é $y_i$. A log-verossimilhança dos parâmetros $\beta$ é, portanto:

$$\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 uma única amostra de dados

Para simplificar a notação, vamos considerar primeiro o caso de um único ponto de dados, $N=1$; em seguida, vamos estender para o caso geral por aditividade.

Gradiente

Temos:

$$\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*}$$

Portanto, pela regra da cadeia:

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

Separadamente, de acordo com "Média e variância da estatística suficiente", temos $A'(\theta) = {\text{Mean}_T}(x^\top \beta)$. Portanto, de acordo com "Afirmação: expressando $h'$ em termos da estatística suficiente", temos:

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

Hessiana

Fazendo a diferenciação uma segunda vez pela regra do produto, obtemos:

$$\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*}$$

Informação de Fisher

De acordo com "Média e variância da estatística suficiente", temos:

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

Portanto:

$$\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 diversas amostras de dados

Agora, estendemos o caso de $N=1$ para o caso geral. Seja $\boldsymbol{\eta} := \mathbf{x} \beta$ o vetor cuja $i$ésima coordenada é a resposta linear da $i$ésima amostra de dados. Seja $\mathbf{T}$ (respectivamente, ${\textbf{Mean}_T}$, respectivamente, ${\textbf{Var}_T}$) a função (vetorizada) difundida que aplica a função de valor escalar $T$ (respectivamente, ${\text{Mean}_T}$, respectivamente, ${\text{Var}_T}$) a cada coordenada. Então, temos:

$$\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*}$$

e

$$\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*}$$

em que as frações denotam a divisão com elementos.

Verificação das fórmulas numericamente

Agora, verificamos numericamente a fórmula acima para o gradiente da log-verossimilhança usando tf.gradients, e verificamos a fórmula para a informação de Fisher com uma estimativa de Monte Carlo usando 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]]

Referências

[1]: Guo-Xun Yuan, Chia-Hua Ho e Chih-Jen Lin. An Improved GLMNET for L1-regularized Logistic Regression (GLMNET aprimorada para regressão logística com regularização L1). Journal of Machine Learning Research, 13, 2012. http://www.jmlr.org/papers/volume13/yuan12a/yuan12a.pdf

[2]: skd. Derivation of Soft Thresholding Operator (Derivação de operador de soft threshold). 2018. https://math.stackexchange.com/q/511106

[3]: Contribuidores da Wikipedia. Proximal gradient methods for learning (Métodos do gradiente proximal para aprendizado). Wikipedia, The Free Encyclopedia, 2018. https://en.wikipedia.org/wiki/Proximal_gradient_methods_for_learning

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