Mô hình tuyến tính tổng quát

Xem trên TensorFlow.org Chạy trong Google Colab Xem nguồn trên GitHubTải xuống sổ ghi chép

Trong sổ tay này, chúng tôi giới thiệu Mô hình tuyến tính tổng quát thông qua một ví dụ đã làm việc. Chúng tôi giải quyết ví dụ này theo hai cách khác nhau bằng cách sử dụng hai thuật toán để điều chỉnh GLM một cách hiệu quả trong Xác suất TensorFlow: chấm điểm Fisher cho dữ liệu dày đặc và giảm độ dốc gần theo tọa độ cho dữ liệu thưa thớt. Chúng tôi so sánh các hệ số được trang bị với hệ số đúng và, trong trường hợp của coordinatewise gần gradient descent, với sản lượng tương tự R của glmnet thuật toán. Cuối cùng, chúng tôi cung cấp thêm các chi tiết toán học và dẫn xuất của một số thuộc tính chính của GLM.

Lý lịch

Một mô hình tuyến tính tổng quát (GLM) là một mô hình tuyến tính (\(\eta = x^\top \beta\)) bọc trong một chuyển đổi (chức năng liên kết) và được trang bị với một phân phối phản hồi từ một gia đình theo cấp số nhân. Việc lựa chọn chức năng liên kết và phân phối phản hồi là rất linh hoạt, điều này cho phép các GLM có tính biểu hiện cao. Chi tiết đầy đủ, bao gồm trình bày tuần tự tất cả các định nghĩa và kết quả xây dựng lên GLM ở dạng ký hiệu rõ ràng, được tìm thấy trong "Nguồn gốc của Sự kiện GLM" bên dưới. Chúng tôi tóm tắt:

Trong một GLM, một phân phối dự đoán cho biến phản ứng \(Y\) được liên kết với một vector các nhân tố ảnh quan sát \(x\). Bản phân phối có dạng:

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

Dưới đây \(\beta\) là các thông số ( "trọng lượng"), \(\phi\) một hyperparameter đại diện phân tán ( "sai"), và \(m\), \(h\), \(T\), \(A\) được đặc trưng bởi các mô hình người dùng chỉ định gia đình.

Giá trị trung bình của \(Y\) phụ thuộc vào \(x\) bởi thành phần của tuyến tính phản ứng \(\eta\) và (ngược lại) chức năng liên kết, ví dụ:

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

nơi \(g\) là cái gọi là chức năng liên kết. Trong TFP lựa chọn chức năng liên kết và mô hình gia đình đang cùng nhau specifed bởi một tfp.glm.ExponentialFamily lớp con. Những ví dụ bao gồm:

TFP thích đặt tên cho gia đình mô hình theo sự phân bố trên Y chứ không phải là chức năng liên kết từ tfp.Distribution s đã là công dân hạng nhất. Nếu tfp.glm.ExponentialFamily lớp con tên chứa một từ thứ hai, điều này cho thấy một chức năng liên kết không kinh điển .

GLM có một số đặc tính đáng chú ý cho phép thực hiện hiệu quả công cụ ước tính khả năng xảy ra tối đa. Đứng đầu trong số các đặc tính này là công thức đơn giản cho gradient của loga \(\ell\), và cho ma trận thông tin Fisher, đó là giá trị kỳ vọng của Hessian của âm loga dưới một tái lấy mẫu của phản ứng dưới cùng một dự báo. I E:

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

nơi \(\mathbf{x}\) là ma trận mà \(i\)thứ hàng là vector dự đoán cho \(i\)thứ mẫu dữ liệu, và \(\mathbf{y}\) là vector có \(i\)thứ phối hợp là phản ứng quan sát cho \(i\)thứ mẫu dữ liệu . Ở đây (lỏng lẻo nói), \({\text{Mean}_T}(\eta) := \mathbb{E}[T(Y)\,|\,\eta]\) và \({\text{Var}_T}(\eta) := \text{Var}[T(Y)\,|\,\eta]\), và in đậm biểu thị vector của các chức năng này. Bạn có thể tìm thấy chi tiết đầy đủ về những phân phối mà những kỳ vọng và sự khác biệt này đã qua trong "Nguồn gốc của Sự kiện GLM" bên dưới.

Một ví dụ

Trong phần này chúng tôi mô tả ngắn gọn và giới thiệu hai built-in GLM phù hợp thuật toán trong TensorFlow Xác suất: Fisher điểm ( tfp.glm.fit ) và coordinatewise gần gradient descent ( tfp.glm.fit_sparse ).

Tập dữ liệu tổng hợp

Hãy giả vờ tải một số tập dữ liệu đào tạo.

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

Lưu ý: Kết nối với thời gian chạy cục bộ.

Trong sổ tay này, chúng tôi chia sẻ dữ liệu giữa các nhân Python và R bằng cách sử dụng các tệp cục bộ. Để kích hoạt tính năng chia sẻ này, vui lòng sử dụng thời gian chạy trên cùng một máy mà bạn có quyền đọc và ghi các tệp cục bộ.

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

Không có quy định L1

Chức năng tfp.glm.fit cụ Fisher ghi bàn, trong đó có như một số các đối số của nó:

  • model_matrix = \(\mathbf{x}\)
  • response = \(\mathbf{y}\)
  • model = callable đó, được đưa ra tranh luận \(\boldsymbol{\eta}\), trả về $ triple \ left ({\ textbf {Mean} _T} (\ boldsymbol {\ eta}) {\ textbf {Var} _T} (\ boldsymbol {\ eta} ), {\ textbf {Mean} _T} '(\ boldsymbol {\ eta}) \ right) $.

Chúng tôi khuyến cáo model là một thể hiện của các tfp.glm.ExponentialFamily lớp. Có một số triển khai được tạo sẵn có sẵn, vì vậy đối với hầu hết các GLM phổ biến, không cần mã tùy chỉnh.

@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

Chi tiết toán học

Tính điểm Fisher là một sửa đổi của phương pháp Newton để tìm ước tính khả năng xảy ra tối đa

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

Phương pháp của Vanilla Newton, tìm kiếm các số không của gradient của khả năng xảy ra nhật ký, sẽ tuân theo quy tắc cập nhật

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

nơi \(\alpha \in (0, 1]\) là một tỷ lệ học dùng để điều khiển kích thước bước.

Trong tính điểm Fisher, chúng tôi thay thế Hessian bằng ma trận thông tin Fisher phủ định:

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

[Lưu ý rằng đây \(\mathbf{Y} = (Y_i)_{i=1}^{n}\) là ngẫu nhiên, trong khi \(\mathbf{y}\) vẫn là vector của các phản ứng quan sát được.]

Bằng các công thức trong "Phù hợp các tham số GLM với dữ liệu" bên dưới, điều này đơn giản hóa thành

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

Với sự điều tiết L1

tfp.glm.fit_sparse cụ một GLM fitter hơn phù hợp với các bộ dữ liệu thưa thớt, dựa trên các thuật toán trong Yuan, Ho và Lin 2012 . Các tính năng của nó bao gồm:

  • Chính quy hóa L1
  • Không có nghịch đảo ma trận
  • Vài đánh giá về gradient và Hessian.

Đầu tiên chúng tôi trình bày một ví dụ về cách sử dụng mã. Chi tiết về các thuật toán được nói rõ hơn trong "Chi tiết thuật toán cho tfp.glm.fit_sparse " dưới đây.

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:

Lưu ý rằng các hệ số đã học có cùng dạng thưa thớt với các hệ số thực.

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

So với R của glmnet

Chúng tôi so sánh đầu ra của coordinatewise gần gradient descent với R của glmnet , trong đó sử dụng một thuật toán tương tự.

LƯU Ý: Để thực hiện phần này, bạn phải chuyển sang thời gian chạy 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)

So sánh hệ số R, TFP và true (Lưu ý: quay lại hạt nhân 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

Thông tin chi tiết thuật toán cho tfp.glm.fit_sparse

Chúng tôi trình bày thuật toán dưới dạng một chuỗi ba sửa đổi đối với phương pháp của Newton. Trong mỗi một, nguyên tắc cập nhật \(\beta\) được dựa trên một vector \(s\) và một ma trận \(H\) mà xấp xỉ gradient và Hessian của loga. Trong bước \(t\), chúng tôi chọn một phối hợp \(j^{(t)}\) để thay đổi, và chúng tôi cập nhật \(\beta\) theo nguyên tắc cập nhật:

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

Bản cập nhật này là một Newton-như bước với học tỷ lệ \(\alpha\). Ngoại trừ các mảnh thức (L1 quy tắc), những sửa đổi dưới đây chỉ khác nhau ở cách họ cập nhật \(s\) và \(H\).

Điểm xuất phát: Phương pháp tọa độ Newton's

Trong phương pháp coordinatewise Newton, chúng tôi đặt \(s\) và \(H\) để gradient đúng và Hessian của loga:

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

Ít đánh giá hơn về gradient và Hessian

Gradient và Hessian của khả năng log thường rất tốn kém để tính toán, vì vậy việc tính gần đúng chúng thường rất đáng giá. Chúng ta có thể làm như sau:

  • Thông thường, tính gần đúng Hessian là hằng số cục bộ và xấp xỉ gradient theo bậc đầu tiên bằng cách sử dụng Hessian (gần đúng):

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

  • Đôi khi thực hiện một "vanilla" cập nhật bước như trên, thiết \(s^{(t+1)}\) để gradient chính xác và \(H^{(t+1)}\) đến Hessian chính xác của loga, đánh giá ở \(\beta^{(t+1)}\).

Thay thế thông tin tiêu cực của Fisher cho Hessian

Để tiếp tục giảm chi phí của các bước cập nhật vani, chúng ta có thể thiết lập \(H\) vào ma trận thông tin tiêu cực Fisher (hiệu quả tính toán bằng cách sử dụng công thức trong "Lắp GLM thông số để dữ liệu" bên dưới) chứ không phải là chính xác 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*} \]

Điều hòa L1 thông qua chuyển tiếp dốc gần

Để kết hợp chính quy L1, chúng tôi thay thế quy tắc cập nhật

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

với quy tắc cập nhật chung hơn

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

nơi \(r_{\text{L1} } > 0\) là một hằng số đã cung cấp (L1 chính quy hệ số) và \(\text{SoftThreshold}\) là các nhà điều hành ngưỡng mềm, được xác định bởi

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

Quy tắc cập nhật này có hai thuộc tính truyền cảm hứng sau đây, chúng tôi giải thích bên dưới:

  1. Trong giới hạn trường hợp \(r_{\text{L1} } \to 0\) (tức là, không có quy tắc L1), quy tắc bản cập nhật này là giống với nguyên tắc cập nhật ban đầu.

  2. Quy tắc cập nhật này có thể được hiểu là áp dụng toán tử lân cận có điểm cố định là giải pháp cho vấn đề tối thiểu hóa điều chỉnh 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). $$

Trường hợp thoái hóa \(r_{\text{L1} } = 0\) phục hồi nguyên tắc cập nhật ban đầu

Để xem (1), lưu ý rằng nếu \(r_{\text{L1} } = 0\) sau đó \(\gamma^{(t)} = 0\), do đó

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

Kể từ đây

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

Toán tử lân cận có điểm cố định là MLE chính quy

Để xem (2), nốt nhạc đầu tiên (xem Wikipedia ) mà đối với bất kỳ \(\gamma > 0\), các quy tắc cập nhật

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

đáp ứng (2), nơi \(\text{prox}\) là các nhà điều hành tiệm cận (xem Yu , nơi điều hành này được ký hiệu \(\mathsf{P}\)). Phía bên tay phải của phương trình trên được tính ở đây :

$$

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

Đặc biệt, thiết\(\gamma = \gamma^{(t)} = -\frac{\alpha\, r_{\text{L1} } }{\left(H^{(t)}\right)_{j^{(t)}, j^{(t)} } }\)(lưu ý rằng \(\gamma^{(t)} > 0\) miễn là loga tiêu cực là lồi), chúng tôi có được sự cai trị cập nhật

$$

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

Sau đó chúng tôi thay thế gradient chính xác $ \ left (\ nabla \ beta \, \ ell (\ beta \,; \, \ mathbf {x}, \ mathbf {y}) \ right) {\ beta = \ beta ^ {( t)}} $ với xấp xỉ của nó \(s^{(t)}\), thu thập

\ begin {class} \ left (\ beta {\ text {chính xác prox}, \ gamma ^ {(t)}} ^ {(t + 1)} \ right) {j ^ {(t)}} & \ xấp xỉ \ 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)} \bên phải). \ end {class}

Kể từ đây

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

Nguồn gốc của sự kiện GLM

Trong phần này, chúng tôi trình bày đầy đủ chi tiết và lấy kết quả về GLM được sử dụng trong các phần trước. Sau đó, chúng tôi sử dụng TensorFlow của gradients để số lượng xác minh các công thức có nguồn gốc cho gradient của loga và thông tin Fisher.

Điểm và thông tin về Fisher

Hãy xem xét một gia đình phân bố xác suất tham số của véc tơ tham số \(\theta\), có mật độ xác xuất \(\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T} }\). Điểm số của một kết quả \(y\) tại tham số vector \(\theta_0\) được định nghĩa là độ chênh lệch của log likelihood của \(y\) (đánh giá ở \(\theta_0\)), có nghĩa là,

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

Yêu cầu: Kỳ vọng về điểm số là 0

Trong các điều kiện bình thường nhẹ (cho phép chúng ta vượt qua sự khác biệt theo tích phân),

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

Bằng chứng

Chúng ta có

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

trong đó chúng tôi đã sử dụng: (1) quy tắc chuỗi để phân biệt, (2) định nghĩa kỳ vọng, (3) vượt qua phân biệt dưới dấu tích phân (sử dụng các điều kiện chính quy), (4) tích phân của mật độ xác suất là 1.

Tuyên bố (thông tin Fisher): Phương sai của điểm bằng Hessian dự kiến ​​phủ định của khả năng nhật ký

Trong các điều kiện bình thường nhẹ (cho phép chúng ta vượt qua sự khác biệt theo tích phân),

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

nơi \(\nabla_\theta^2 F\) biểu thị ma trận Hessian, mà \((i, j)\) entry là \(\frac{\partial^2 F}{\partial \theta_i \partial \theta_j}\).

Phía bên tay trái của phương trình này được gọi là thông tin Fisher của gia đình \(\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T} }\) tại vector tham số \(\theta_0\).

Bằng chứng xác nhận quyền sở hữu

Chúng ta có

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

trong đó chúng tôi đã sử dụng (1) quy tắc chuỗi để phân biệt, (2) quy tắc thương cho sự khác biệt, (3) lại quy tắc chuỗi, ngược lại.

Để hoàn thành bằng chứng, chỉ cần chứng minh rằng

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

Để làm điều đó, chúng ta vượt qua sự phân biệt dưới dấu tích phân hai lần:

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

Bổ đề về đạo hàm của hàm phân hoạch log

Nếu \(a\), \(b\) và \(c\) là hàm vô hướng có giá trị, \(c\) hai lần khả vi, như vậy mà gia đình các bản phân phối \(\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T} }\) xác định bởi

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

đáp ứng các điều kiện đều đặn nhẹ mà cho phép đi qua sự khác biệt liên quan đến \(\theta\) dưới một thể thiếu đối với các \(y\), sau đó

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

và

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

(Đây \('\) biểu thị sự khác biệt, vì vậy \(c'\) và \(c''\) là đạo hàm thứ nhất và thứ hai của \(c\).)

Bằng chứng

Đối với gia đình này phân phối, chúng tôi có \(\text{score}(y, \theta_0) = b(y) - c'(\theta_0)\). Phương trình đầu tiên sau đó sau từ thực tế là \(\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ \text{score}(y, \theta_0) \right] = 0\). Tiếp theo, chúng tôi có

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

Gia đình lũy thừa phân tán

A (vô hướng) overdispersed gia đình mũ là một gia đình của các bản phân phối có mật độ mang hình thức

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

nơi \(m\) và \(T\) là được biết đến chức năng vô hướng có giá trị, và \(\theta\) và \(\phi\) là các thông số vô hướng.

[Lưu ý rằng \(A\) được overdetermined: đối với bất kỳ \(\phi_0\), hàm \(A\) là hoàn toàn xác định bởi các hạn chế mà\(\int p_{\text{OEF}(m, T)}(y\ |\ \theta, \phi=\phi_0)\, dy = 1\)cho tất cả \(\theta\). Các \(A\)'s sản xuất bởi các giá trị khác nhau của \(\phi_0\) tất cả phải giống nhau, mà đặt một hạn chế về chức năng \(m\) và \(T\).]

Trung bình và phương sai của thống kê đủ

Trong các điều kiện tương tự như "Bổ đề về đạo hàm của hàm phân hoạch log", chúng ta có

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

\right]

A'(\theta) $$

và

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

\right]

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

Bằng chứng

Theo "Bổ đề về đạo hàm của hàm phân hoạch log", chúng ta có

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

\right]

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

và

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

\right]

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

Kết quả sau đó sau từ thực tế là kỳ vọng là tuyến tính (\(\mathbb{E}[aX] = a\mathbb{E}[X]\)) và phương sai là độ-2 đồng nhất (\(\text{Var}[aX] = a^2 \,\text{Var}[X]\)).

Mô hình tuyến tính tổng quát

Trong một mô hình tuyến tính tổng quát, một phân phối dự đoán cho biến phản ứng \(Y\) được liên kết với một vector các nhân tố ảnh quan sát \(x\). Sự phân bố là thành viên của một gia đình mũ overdispersed, và tham số \(\theta\) được thay thế bằng \(h(\eta)\) nơi \(h\) là một chức năng được biết đến, \(\eta := x^\top \beta\) là cái gọi là phản ứng tuyến tính, và \(\beta\) là một vector của các tham số (hệ số hồi quy) được học. Nói chung các tham số phân tán \(\phi\) có thể được học quá, nhưng trong thiết lập của chúng tôi, chúng tôi sẽ đối xử với \(\phi\) như được biết đến. Vì vậy, thiết lập của chúng tôi là

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

nơi cấu trúc mô hình được đặc trưng bởi sự phân bố \(p_{\text{OEF}(m, T)}\) và chức năng \(h\) mà chuyển đổi đáp ứng tuyến tính với các thông số.

Theo truyền thống, các ánh xạ từ tuyến tính phản ứng \(\eta\) để bình \(\mu := \mathbb{E}_{Y \sim p_{\text{OEF}(m, T)}(\cdot\, |\, \theta = h(\eta), \phi)}\left[ Y\right]\) được ký hiệu

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

Lập bản đồ này là cần thiết để có một-một, và nghịch đảo của nó, \(g\), được gọi là chức năng liên kết cho GLM này. Thông thường, người ta mô tả một GLM bằng cách đặt tên cho hàm liên kết và họ phân phối của nó - ví dụ: "GLM với phân phối Bernoulli và hàm liên kết logit" (còn được gọi là mô hình hồi quy logistic). Để mô tả đầy đủ các GLM, hàm \(h\) cũng phải được xác định. Nếu \(h\) là bản sắc, sau đó \(g\) được cho là chức năng liên kết kinh điển.

Khẳng định: Bày tỏ \(h'\) về số liệu thống kê đầy đủ

Định nghĩa

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

và

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

Sau đó chúng tôi có

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

Bằng chứng

Theo "Trung bình và phương sai của thống kê đủ", chúng tôi có

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

Phân biệt với quy tắc chuỗi, chúng tôi có được

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

và bằng "Trung bình và phương sai của thống kê đủ",

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

Kết luận sau đây.

Phù hợp các thông số GLM với dữ liệu

Các tính chất có nguồn gốc ở trên cho vay mình rất tốt để phù hợp các thông số GLM \(\beta\) đến một tập dữ liệu. Các phương pháp Quasi-Newton chẳng hạn như tính điểm Fisher dựa vào độ dốc của khả năng xảy ra nhật ký và thông tin Fisher, mà hiện nay chúng tôi hiển thị có thể được tính toán đặc biệt hiệu quả cho GLM.

Giả sử chúng ta có vectơ yếu tố dự báo quan sát \(x_i\) và phản ứng vô hướng liên \(y_i\). Ở dạng ma trận, chúng tôi sẽ nói rằng chúng ta có dự đoán quan sát \(\mathbf{x}\) và phản ứng \(\mathbf{y}\), nơi \(\mathbf{x}\) là ma trận mà \(i\)thứ hàng là \(x_i^\top\) và \(\mathbf{y}\) là vector có \(i\)thứ phần tử là \(y_i\). Khả năng ghi các thông số \(\beta\) là sau đó

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

Đối với một mẫu dữ liệu duy nhất

Để đơn giản hóa các ký hiệu, trước tiên hãy xem xét trường hợp của một điểm dữ liệu duy nhất, \(N=1\); thì chúng ta sẽ mở rộng đến trường hợp chung bằng tính cộng.

Dốc

Chúng ta có

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

Do đó theo quy tắc chuỗi,

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

Riêng biệt, bởi "Mean và phương sai của các số liệu thống kê đầy đủ," chúng tôi có \(A'(\theta) = {\text{Mean}_T}(x^\top \beta)\). Do đó, bằng cách "yêu cầu bồi thường: Bày tỏ \(h'\) về số liệu thống kê đầy đủ," chúng tôi có

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

Hessian

Phân biệt lần thứ hai, theo quy tắc sản phẩm, chúng tôi có được

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

Thông tin về Fisher

Theo "Trung bình và phương sai của thống kê đủ", chúng tôi có

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

Kể từ đây

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

Đối với nhiều mẫu dữ liệu

Bây giờ chúng ta mở rộng \(N=1\) vụ án cho trường hợp tổng quát. Hãy \(\boldsymbol{\eta} := \mathbf{x} \beta\) biểu thị vector mà \(i\)thứ phối hợp là phản ứng tuyến tính từ \(i\)thứ mẫu dữ liệu. Hãy \(\mathbf{T}\) (resp. \({\textbf{Mean}_T}\), resp. \({\textbf{Var}_T}\)) biểu thị phát sóng chức năng (vector) mà áp dụng vô hướng có giá trị chức năng \(T\) (resp. \({\text{Mean}_T}\), resp. \({\text{Var}_T}\)) cho mỗi phối hợp . Sau đó chúng tôi có

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

và

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

trong đó các phân số biểu thị phép chia theo nguyên tố.

Xác minh công thức bằng số

Bây giờ chúng ta xác minh công thức trên cho gradient của log likelihood số lượng sử dụng tf.gradients , và xác minh công thức để biết thông tin Fisher với một ước tính Monte Carlo sử dụng 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]]

Người giới thiệu

[1]: Guo-Xun Yuan, Chia-Hua Ho và Chih-Jen Lin. GLMNET được cải tiến cho hồi quy logistic chính quy L1. Tạp chí Nghiên cứu Machine Learning, 13, 2012. http://www.jmlr.org/papers/volume13/yuan12a/yuan12a.pdf

[2]: skd. Khởi nguồn của toán tử ngưỡng mềm. 2018. https://math.stackexchange.com/q/511106

[3]: Người đóng góp Wikipedia. Phương pháp học tập gần bằng gradient. Wikipedia, bách khoa toàn thư miễn phí, năm 2018. https://en.wikipedia.org/wiki/Proximal_gradient_methods_for_learning

[4]: Yao-Liang Yu. Toán tử lân cận. https://www.cs.cmu.edu/~suvrit/teach/yaoliang_proximity.pdf