Посмотреть на TensorFlow.org

Запускаем в Google Colab

Посмотреть исходный код на GitHub

Скачать блокнот

В этой записной книжке мы представляем обобщенные линейные модели на рабочем примере. Мы решаем этот пример двумя разными способами, используя два алгоритма для эффективной подгонки GLM в TensorFlow Probability: оценка Фишера для плотных данных и покоординатный проксимальный градиентный спуск для разреженных данных. Мы сравниваем подогнанные коэффициенты для истинных коэффициентов и, в случае покоординатного проксимальные градиентный спуск, к выходу аналогичного АиРа glmnet алгоритма. Наконец, мы предоставляем дополнительные математические детали и выводы нескольких ключевых свойств GLM.

Задний план

Обобщенная линейная модель (GLM) представляет собой линейную модель (\(\eta = x^\top \beta\)) , завернутые в трансформации (функции связи) и оборудованный с распределением отклика от экспоненциального семейства. Выбор функции связи и распределения ответов очень гибкий, что придает большую выразительность GLM. Полную информацию, включая последовательное представление всех определений и результатов, составляющих GLM, в однозначной нотации, можно найти в разделе «Вывод фактов GLM» ниже. Резюмируем:

В GLM, предиктивные распределения для переменного отклика \(Y\) связанно с вектором наблюдаемых предикторами \(x\). Распределение имеет вид:

\[ \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*} \]

Здесь \(\beta\) являются параметры ( "весы"), \(\phi\) гиперпараметра , представляющей дисперсия ( "дисперсия"), и \(m\), \(h\), \(T\), \(A\) характеризуется указанным пользователем модели семья.

Среднее \(Y\) зависит от \(x\) композицией линейного отклика \(\eta\) и (обратной) линии связи функции, а именно:

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

где \(g\) является так называемой функцией связи. В TFP выбор функции связи и модели семьи совместно с помощью предписанных в технических заданиях tfp.glm.ExponentialFamily подкласса. Примеры включают:

• tfp.glm.Normal , также известная как «линейная регрессия»
• tfp.glm.Bernoulli , также известный как «логистическая регрессия»
• tfp.glm.Poisson , ака «Пуассон регрессия»
• tfp.glm.BernoulliNormalCDF , также известный как «пробит регрессии».

TFP предпочитает называть модели семьи в соответствии с распределением по Y , а не функцию связи , так как tfp.Distribution s уже гражданами первого класса. Если в tfp.glm.ExponentialFamily подкласса имя содержит второе слово, то это указывает на неканоничной функции связи .

GLM обладают несколькими замечательными свойствами, которые позволяют эффективно реализовать оценку максимального правдоподобия. Главная из этих свойств являются простыми формулами для градиента логарифмического правдоподобия \(\ell\), а также для информационной матрицы Фишера, который является ожидаемым значением гессиана отрицательного логарифмического правдоподобия при повторной выборке ответа в соответствии с те же предикторы. Т.е.:

\[ \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*} \]

где \(\mathbf{x}\) является матрицей, \(i\)й строкой является прогностическим вектором для \(i\)- й выборки данных, а также \(\mathbf{y}\) является вектором, \(i\)й координаты является наблюдаемым откликом для \(i\)- й выборки данных . Здесь (грубо говоря), \({\text{Mean}_T}(\eta) := \mathbb{E}[T(Y)\,|\,\eta]\) и \({\text{Var}_T}(\eta) := \text{Var}[T(Y)\,|\,\eta]\)и полужирный обозначает векторизации этих функций. Полную информацию о распределении этих ожиданий и отклонений можно найти в разделе «Получение фактов GLM» ниже.

Пример

В этом разделе мы кратко опишем и продемонстрировать два встроенных GLM установки алгоритмов в TensorFlow Вероятность: Fisher скоринг ( tfp.glm.fit ) и покоординатное проксимального градиент снижения ( tfp.glm.fit_sparse ).

Синтетический набор данных

Давайте представим, что загружаем некоторый обучающий набор данных.

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

Примечание. Подключитесь к локальной среде выполнения.

В этой записной книжке мы обмениваемся данными между ядрами Python и R с помощью локальных файлов. Чтобы включить этот общий доступ, используйте среду выполнения на том же компьютере, на котором у вас есть разрешение на чтение и запись локальных файлов.

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

Без регуляризации L1

Функция tfp.glm.fit реализует Fisher скоринга, который принимает как некоторые из его аргументов:

• model_matrix = \(\mathbf{x}\)
• response = \(\mathbf{y}\)
• model = вызываемым который, учитывая аргумент \(\boldsymbol{\eta}\), возвращает тройной $ \ влево ({\ textbf {Среднее} _t} (\ boldsymbol {\ ETA}), {\ textbf {вар} _t} (\ boldsymbol {\ ETA} ), {\ textbf {Mean} _T} '(\ boldsymbol {\ eta}) \ right) $.

Мы рекомендуем model быть экземпляром tfp.glm.ExponentialFamily класса. Доступно несколько готовых реализаций, поэтому для большинства распространенных GLM специальный код не требуется.

@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

Математические детали

Оценка Фишера - это модификация метода Ньютона для нахождения оценки максимального правдоподобия.

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

Метод Ванили Ньютона, ищущий нули градиента логарифмической вероятности, будет следовать правилу обновления

\[ \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)} } \]

где \(\alpha \in (0, 1]\) является скорость обучения используется для управления размером шага.

В системе оценок Фишера мы заменяем гессиан на отрицательную информационную матрицу Фишера:

\[ \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*} \]

[Следует отметить , что здесь \(\mathbf{Y} = (Y_i)_{i=1}^{n}\) является случайным, тогда как \(\mathbf{y}\) по - прежнему вектор наблюдаемых реакций.]

По формулам, приведенным ниже в разделе «Подбор параметров GLM к данным», это упрощается до

\[ \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*} \]

С регуляризацией L1

tfp.glm.fit_sparse реализует GLM слесарем более подходит для разреженных наборов данных на основе алгоритма в Юань, Хо и Лин 2012 . Его особенности включают в себя:

• L1 регуляризация
• Нет инверсии матриц
• Немногочисленные оценки градиента и гессиана.

Сначала мы представляем пример использования кода. Детали алгоритма доработаны в «Алгоритм Детали для tfp.glm.fit_sparse » ниже.

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:

Обратите внимание, что изученные коэффициенты имеют тот же образец разреженности, что и истинные коэффициенты.

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

Сравните с R - х glmnet

Сравним выход покоординатное проксимальный градиентного спуска к тому , что из АиР glmnet , который использует подобный алгоритм.

ПРИМЕЧАНИЕ. Для выполнения этого раздела необходимо переключиться в среду выполнения R colab.

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)

Сравните R, TFP и истинные коэффициенты (Примечание: назад к ядру 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

Алгоритм Детали для tfp.glm.fit_sparse

Мы представляем алгоритм как последовательность трех модификаций метода Ньютона. В каждом из них, правило обновления для \(\beta\) основан на векторном \(s\) и матрицы \(H\) , который аппроксимирует градиент и гессианом логарифмического правдоподобия. На этапе \(t\), мы выбираем координату \(j^{(t)}\) изменения, и мы обновляем \(\beta\) согласно правилу обновления:

\[ \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*} \]

Это обновление является Newton-подобный шаг с изучением скорости \(\alpha\). Для окончательной части (L1) , за исключением регуляризации, приведенные ниже модификации различаются только в том , как они обновляются \(s\) и \(H\).

Отправная точка: координатный метод Ньютона.

В методе покоординатного Ньютона, мы устанавливаем \(s\) и \(H\) к истинному градиента и гессиана журнала правдоподобия:

\[ \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*} \]

Меньше оценок градиента и гессиана

Градиент и гессиан логарифма правдоподобия часто являются дорогостоящими для вычисления, поэтому часто имеет смысл их аппроксимировать. Сделать это можно следующим образом:

• Обычно аппроксимируют гессиан как локально постоянный и аппроксимируют градиент до первого порядка, используя (приблизительный) гессиан:

\[ \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*} \]

• Иногда выполнить «ванильный» шаг обновления , как указаны выше, установку \(s^{(t+1)}\) к точному градиента и \(H^{(t+1)}\) к точному гессиану логарифмического правдоподобия, измеренному при \(\beta^{(t+1)}\).

Заменить гессенскую отрицательную информацию Фишера

Для дальнейшего снижения стоимости этапов обновления ванили, мы можем установить \(H\) к отрицательному информационной матрице Фишера (эффективно вычислимым с использованием формул в «Установке параметров GLM к данным» ниже) , а не точной Hessian:

\[ \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*} \]

L1 регуляризация через проксимальный градиентный спуск

Чтобы включить регуляризацию L1, мы заменяем правило обновления

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

с более общим правилом обновления

\[ \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*} \]

где \(r_{\text{L1} } > 0\) является постоянным Прилагаемые (коэффициент регуляризации L1) и \(\text{SoftThreshold}\) является мягким оператором порогового значения, определяемого

\[ \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} \]

Это правило обновления имеет следующие два вдохновляющих свойства, которые мы объясним ниже:

1. В предельном случае \(r_{\text{L1} } \to 0\) (т.е. без L1 регуляризации), это правило обновление идентично исходному правилу обновления.

2. Это правило обновления можно интерпретировать как применение оператора близости, фиксированная точка которого является решением 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). $$

Вырожденный случай \(r_{\text{L1} } = 0\) восстанавливает исходное правило обновления

Чтобы увидеть (1), заметим , что если \(r_{\text{L1} } = 0\) затем \(\gamma^{(t)} = 0\), следовательно ,

\[ \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*} \]

Следовательно

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

Оператор близости, неподвижной точкой которого является регуляризованный MLE

Чтобы увидеть (2), первую ноту (см Википедии ) , что для любого \(\gamma > 0\), правило обновления

\[ \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) \]

удовлетворяет условию (2), где \(\text{prox}\) является оператором близости (см Ю , где этот оператор обозначается \(\mathsf{P}\)). Правая часть приведенного выше уравнения вычисляется здесь :

$$

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

В частности, установка\(\gamma = \gamma^{(t)} = -\frac{\alpha\, r_{\text{L1} } }{\left(H^{(t)}\right)_{j^{(t)}, j^{(t)} } }\)(обратите внимание , что \(\gamma^{(t)} > 0\) до тех пор , как отрицательный логарифм правдоподобия выпукло), получаем правило обновления

$$

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

Затем заменим точный градиент $ \ влево (\ набла \ бета \, \ ell_p (\ бета \; \, \ mathbf {X}, \ mathbf {у}) \ справа) {\ бета = \ бета ^ {( т)}} $ с его приближением \(s^{(t)}\), получение

\ начать {Выровнять} \ влево (\ бета {\ текст {точным Prox}, \ Gamma ^ {(т)}} ^ {(т + 1)} \ справа) {J ^ {(т)}} & \ примерно \ текст {SoftThreshold} \ влево (\ бета ^ {(т)} {J ^ {(т)}} - \ альфа \ гидроразрыва {\ влево (s ^ {(т)} \ справа) {J ^ {( т)}}} {\ влево (Н ^ {(т)} \ справа) {J ^ {(т)}, J ^ {(т)}}}, \ \ гамма ^ {(т)} \ справа) \ & = \ текст {SoftThreshold} \ влево (\ бета ^ {(т)} {J ^ {(т)}} - \ альфа \, и ^ {(т)}, \ \ гамма ^ {(т)} \правильно). \ конец {} Align

Следовательно

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

Вывод фактов GLM

В этом разделе мы подробно излагаем и выводим результаты о GLM, которые использовались в предыдущих разделах. Затем мы используем TensorFlow в gradients численно проверить полученные формулы для градиента логарифмической вероятности и информации Фишера.

Оценка и информация Фишера

Рассмотрит семейство вероятностных распределений параметризованных вектор параметров \(\theta\), имеющей плотность вероятности \(\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T} }\). Оценка итогового \(y\) в параметре вектора \(\theta_0\) определяется как градиент журнала вероятности \(y\) (оценивали при \(\theta_0\)), то есть,

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

Утверждение: ожидание результата равно нулю.

В условиях мягкой регулярности (позволяющей пройти дифференцирование под интегралом)

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

Доказательство

У нас есть

\[ \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*} \]

где мы использовали: (1) цепное правило дифференцирования, (2) определение математического ожидания, (3) переходное дифференцирование под знаком интеграла (с использованием условий регулярности), (4) интеграл плотности вероятности равен 1.

Утверждение (информация Фишера): дисперсия оценки равна отрицательному ожидаемому гессиану логарифмической вероятности.

В условиях мягкой регулярности (позволяющей пройти дифференцирование под интегралом)

$$ \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=\theta0)}\left[ \left(\nabla\theta^2 \log p(Y | \theta)\right)_{\theta=\theta_0} \right] $$

где \(\nabla_\theta^2 F\) обозначает матрицу Гессе, чьи \((i, j)\) запись \(\frac{\partial^2 F}{\partial \theta_i \partial \theta_j}\).

Левая часть этого уравнения называется информацией Фишера семейство \(\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T} }\) при векторе параметров \(\theta_0\).

Доказательство претензии

У нас есть

\[ \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*} \]

где мы использовали (1) цепное правило для дифференцирования, (2) правило частного для дифференцирования, (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) } \right] \stackrel{\text{?} }{=} 0. \]

Для этого дважды пройдем дифференцирование под знаком интеграла:

\[ \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*} \]

Лемма о производной логарифмической статистической суммы

Если \(a\), \(b\) и \(c\) являются скалярные функции, \(c\) дважды дифференцируемые, такие , что семейство распределений \(\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T} }\) , определяемой

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

удовлетворяет мягкие условия регулярности , которые позволяют прохождение дифференцирования по \(\theta\) под интегралом по отношению к \(y\), то

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

и

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

(Здесь \('\) означает дифференцирование, поэтому \(c'\) и \(c''\) является первым и вторым производным \(c\).)

Доказательство

Для этого семейства распределений, мы имеем \(\text{score}(y, \theta_0) = b(y) - c'(\theta_0)\). Первое уравнение следует из того факта , что \(\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ \text{score}(y, \theta_0) \right] = 0\). Далее у нас есть

\[ \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*} \]

Сверхдисперсная экспоненциальная семья

А (скалярный) overdispersed экспоненциального семейства представляет собой семейство распределений плотности которых принимает форму

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

где \(m\) и \(T\) известны скалярные функции и \(\theta\) и \(\phi\) являются скалярными параметрами.

[Следует отметить , что \(A\) переопределена: для любого \(\phi_0\), функция \(A\) полностью определяется тем ограничением , что\(\int p_{\text{OEF}(m, T)}(y\ |\ \theta, \phi=\phi_0)\, dy = 1\)для всех \(\theta\). \(A\)«ы производства различных значений \(\phi_0\) все должны быть такими же, что накладывает ограничение на функции \(m\) и \(T\).]

Среднее и дисперсия достаточной статистики

При тех же условиях, что и «Лемма о производной логарифмической статистической суммы», мы имеем

$$ \mathbb{E}{Y \sim p{\text{OEF}(m, T)}(\cdot | \theta, \phi)} \left[ T(Y)

\right]

A'(\theta) $$

и

$$ \text{Var}{Y \sim p{\text{OEF}(m, T)}(\cdot | \theta, \phi)} \left[ T(Y)

\right]

\phi A''(\theta). $$

Доказательство

По «лемме о производной логарифмической статистической суммы» имеем

$$ \mathbb{E}{Y \sim p{\text{OEF}(m, T)}(\cdot | \theta, \phi)} \left[ \frac{T(Y)}{\phi}

\right]

\frac{A'(\theta)}{\phi} $$

и

$$ \text{Var}{Y \sim p{\text{OEF}(m, T)}(\cdot | \theta, \phi)} \left[ \frac{T(Y)}{\phi}

\right]

\frac{A''(\theta)}{\phi}. $$

Результат следует из того факта , что ожидание является линейным (\(\mathbb{E}[aX] = a\mathbb{E}[X]\)) и дисперсией является степенью 2 однородного (\(\text{Var}[aX] = a^2 \,\text{Var}[X]\)).

Обобщенная линейная модель

В обобщенной линейной модели, предиктивные распределения для переменного отклика \(Y\) связанно с вектором наблюдаемых предикторов \(x\). Распределение является членом overdispersed экспоненциального семейства, а параметр \(\theta\) заменяется \(h(\eta)\) , где \(h\) является известной функцией, \(\eta := x^\top \beta\) является так называемым линейным откликом, и \(\beta\) представляет собой вектор параметров (коэффициентов регрессии), которые необходимо изучить. В общем случае параметр дисперсии \(\phi\) можно научиться тоже, но в нашей установке мы будем рассматривать \(\phi\) , как известно. Итак, наша установка

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

где структура модель характеризуется распределением \(p_{\text{OEF}(m, T)}\) и функции \(h\) , который преобразует линейный отклик на параметры.

Традиционно, отображение из линейного отклика \(\eta\) к среднему \(\mu := \mathbb{E}_{Y \sim p_{\text{OEF}(m, T)}(\cdot\, |\, \theta = h(\eta), \phi)}\left[ Y\right]\) обозначается

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

Это отображение требуется , чтобы один-к-одному, и обратное, \(g\), называется функцией связи для этого GLM. Обычно GLM описывают, называя его функцию связи и семейство распределений - например, «GLM с распределением Бернулли и функцией логита связи» (также известной как модель логистической регрессии). Для того , чтобы полностью охарактеризовать GLM, функция \(h\) также должен быть указан. Если \(h\) является тождественным, то \(g\) называется канонической функцией связи.

Требование: Выражение \(h'\) в терминах достаточной статистики

Определять

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

и

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

Тогда у нас есть

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

Доказательство

Под «средним значением и дисперсией достаточной статистики» мы имеем

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

Дифференцируя цепным правилом, получаем

\[ {\text{Mean}_T}'(\eta) = A''(h(\eta))\, h'(\eta), \]

и "Среднее значение и дисперсия достаточной статистики"

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

Напрашивается вывод.

Подгонка параметров GLM к данным

Свойства полученных выше поддаются очень хорошо подгоночных параметров GLM \(\beta\) к набору данных. Квазиньютоновские методы, такие как оценка Фишера, основываются на градиенте логарифмической вероятности и информации Фишера, которую, как мы теперь показываем, можно вычислить особенно эффективно для GLM.

Предположим , что мы наблюдали векторы прогнозирующих \(x_i\) и связанные с ними скалярные ответов \(y_i\). В матричной форме, мы будем говорить , что мы наблюдали предсказатель \(\mathbf{x}\) и ответ \(\mathbf{y}\), где \(\mathbf{x}\) является матрицей, \(i\)й строкой является \(x_i^\top\) и \(\mathbf{y}\) является вектором, \(i\)- й элемент \(y_i\). Логарифмическое правдоподобие параметров \(\beta\) затем

\[ \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). \]

Для единичной выборки данных

Для упрощения обозначений, давайте сначала рассмотрим случай одной точки данных, \(N=1\); затем мы расширим его на общий случай по аддитивности.

Градиент

У нас есть

\[ \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*} \]

Следовательно, по цепному правилу

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

Отдельно от «среднего и дисперсии достаточной статистики,» мы имеем \(A'(\theta) = {\text{Mean}_T}(x^\top \beta)\). Следовательно, по «претензии: Выражая \(h'\) в терминах достаточной статистики,» мы имеем

\[ \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. \]

Гессен

Дифференцируя второй раз, по правилу произведения получаем

\[ \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*} \]

Информация Fisher

Под «средним значением и дисперсией достаточной статистики» мы имеем

\[ \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. \]

Следовательно

\[ \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*} \]

Для нескольких выборок данных

Продолжим теперь \(N=1\) дело в общем случае. Пусть \(\boldsymbol{\eta} := \mathbf{x} \beta\) обозначим вектор, \(i\)- й координате линейный отклик от \(i\)- й выборки данных. Пусть \(\mathbf{T}\) (соотв. \({\textbf{Mean}_T}\), соответственно \({\textbf{Var}_T}\)) обозначает функцию , которая применяет скалярная функция транслируемое (векторизованную) \(T\) (соответственно \({\text{Mean}_T}\), соответственно \({\text{Var}_T}\)) каждой координате . Тогда у нас есть

\[ \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*} \]

и

\[ \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*} \]

где дроби обозначают поэлементное деление.

Проверка формул численно

Проверим теперь приведенную выше формулу для градиента логарифмического правдоподобия численно с использованием tf.gradients и проверить формулу для информации Фишера с оценкой методом Монте - Карло с использованием 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]]

использованная литература

[1]: Го-Сюнь Юань, Чиа-Хуа Хо и Чжи-Джен Линь. Улучшенный GLMNET для L1-регуляризованной логистической регрессии. Журнал Machine Learning Research, 13, 2012. http://www.jmlr.org/papers/volume13/yuan12a/yuan12a.pdf

[2]: skd. Вывод оператора мягкого порога. 2018. https://math.stackexchange.com/q/511106

[3]: Авторы Википедии. Проксимальные градиентные методы обучения. Википедия, свободная энциклопедия, 2018. https://en.wikipedia.org/wiki/Proximal_gradient_methods_for_learning

[4]: Яо-Лян Юй. Оператор близости. https://www.cs.cmu.edu/~suvrit/teach/yaoliang_proximity.pdf