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.

일반화된 선형 모델

TensorFlow.org에서 보기

Google Colab에서 실행하기

GitHub에서 소그 보기

노트북 다운론드하기

이 노트북에서는 작업 예제를 통해 일반화된 선형 모델을 소개합니다. TensorFlow Probability에서 GLM을 효율적으로 피팅하기 위한 두 가지 알고리즘을 사용하여 이 예제를 두 가지 다른 방법으로 해결합니다. 즉, 밀집 데이터에 대한 Fisher 스코어링과 희소 데이터에 대한 좌표별 근위 경사 하강법입니다. 피팅된 계수를 실제 계수와 비교하고, 좌표별 근위 경사 하강법의 경우 R의 유사한 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$의 평균은 선형 응답 $\eta$와 (역) 연결 함수의 구성에 의해 $x$에 따라 달라집니다. 예:

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

여기서 $g$는 소위 연결 함수입니다. TFP에서 연결 함수와 모델군의 선택은 tfp.glm.ExponentialFamily 하위 클래스에 의해 공동으로 지정됩니다. 예를 들면 다음과 같습니다.

• tfp.glm.Normal, 일명 "선형 회귀"

• tfp.glm.Bernoulli, 일명 "로지스틱 회귀"

• tfp.glm.Poisson, 일명 "푸아송 회귀"

• tfp.glm.BernoulliNormalCDF, 일명 "프로빗 회귀".

tfp.Distribution이 이미 일급 객체이기 때문에 TFP는 연결 함수보다 Y에 대한 분포에 따라 모델군의 이름을 지정하는 것을 선호합니다. tfp.glm.ExponentialFamily 서브 클래스 이름에 두 번째 단어가 포함된 경우, 이는 비표준 연결 함수를 나타냅니다.

GLM에는 최대 가능도 estimator를 효율적으로 구현할 수 있게 해주는 몇 가지 놀라운 속성이 있습니다. 이러한 속성 중 가장 중요한 것은 로그-가능도 $\ell$의 기울기 및 Fisher 정보 행렬에 대한 간단한 공식입니다. 이는 동일한 예측 변수 아래 응답의 재표본 추출에서 음의 로그-가능도에 대한 Hessian 기대값입니다. 즉, 다음과 같습니다.

$$\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$th 행이 $i$th 데이터 샘플에 대한 예측 벡터인 행렬이고 $\mathbf{y}$는 $i$th 좌표가 $i$th 데이터 샘플에 대해 관찰된 응답인 벡터입니다. 여기서(대략적으로 말하면), ${\text{Mean}_T}(\eta) := \mathbb{E}[T(Y),|,\eta]$ 및 ${\text{Var}_T}(\ eta) := \text{Var}[T(Y),|,\eta]$, 및 굵은 글씨체는 이러한 함수의 벡터화를 나타냅니다. 이러한 기대치와 편차가 어떤 분포에 걸쳐져 있는지에 대한 자세한 내용은 아래의 "GLM 사실 도출"에서 확인할 수 있습니다.

예제

이 섹션에서는 TensorFlow Probability의 두 가지 기본 제공 GLM 피팅 알고리즘인 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}$ 인수가 주어지면 삼중 $\left( {\textbf{Mean}_T}(\boldsymbol{\eta}), {\textbf{Var}_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

수학적 세부 사항

Fisher 스코어링은 최대 가능도 추정치를 찾기 위한 Newton의 방법을 변형한 것입니다.

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

로그-가능도의 기울기 0을 찾는 Vanilla Newton의 방법은 업데이트 규칙을 따릅니다.

$$\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]$은 단계 크기를 제어하는 데 사용되는 학습률입니다.

Fisher 스코어링에서는 Hesian을 음의 Fisher 정보 행렬로 바꿉니다.

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

[여기서 $\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는 Yuan, Ho 및 Lin 2012의 알고리즘을 기반으로 희소 데이터세트에 더 적합한 GLM 피터를 구현합니다. 해당 요소에는 다음이 포함됩니다.

• L1 정규화

• 매트릭스 반전 없음

• 기울기와 Hessian에 대한 평가가 거의 없음

먼저 코드의 사용 예를 제시합니다. 알고리즘의 세부 사항은 아래 "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와 비교하기

좌표별 근위 경사 하강법의 출력을 유사한 알고리즘을 사용하는 R의 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에 대한 알고리즘 세부 정보

알고리즘을 Newton의 방법에 대한 세 가지 수정 시퀀스로 제시합니다. 각각에서 $\beta$에 대한 업데이트 규칙은 벡터 $s$와 로그-가능도의 그래디언트 및 Hessian을 근사 계산하는 행렬 $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*}$$

이 업데이트는 학습률이 $\alpha$인 Newton과 유사한 단계입니다. 최종 조각(L1 정규화)을 제외하고 아래 수정 사항은 $s$ 및 $H$를 업데이트하는 방식만 다릅니다.

출발점: 좌표별 Newton 방법

좌표별 Newton 방법에서는 $s$ 및 $H$를 로그-가능도의 실제 그래디언트 및 Hessian으로 설정합니다.

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

그래디언트 및 Hessian에 대한 더 적은 평가

로그-가능도의 그래디언트와 Hessian은 종종 계산하는 데 비용이 많이 들기 때문에 대략적으로 계산하는 것이 좋습니다. 다음과 같이 할 수 있습니다.

• 일반적으로 Hessian을 로컬 상수로 근사하고 (근사) Hessian을 사용하여 그래디언트를 1차로 근사합니다.

$$\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)}$를 로그-가능도의 정확한 Hessian으로 설정합니다. $\beta^{(t+1)}$에서 평가됩니다.

Hessian을 음의 Fisher 정보로 대체

바닐라 업데이트 단계의 비용을 더 줄이기 위해 정확한 Hessian 대신 $H$를 음의 Fisher 정보 행렬(아래 "GLM 매개변수를 데이터에 맞추기"의 공식을 사용하여 효율적으로 계산할 수 있음)로 설정할 수 있습니다.

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

여기서 ParseError: KaTeX parse error: Expected 'EOF', got '&' at position 15: r_{\text{l0}} &̲gt; 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{l0}} \to 0$의 제한적인 경우(즉, L1 정규화 없음), 이 업데이트 규칙은 원래 업데이트 규칙과 동일합니다.

2. $r_{\text{l0}} \to 0$의 제한적인 경우(즉, 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{l0}} = 0$는 원래 업데이트 규칙을 복구합니다.

(1)을 보려면 $r_{\text{l0}} = 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$에 대해 다음 업데이트 규칙(Wikipedia 참조)

$$\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}$로 표기된 Yu 참조). 위 수식의 우변은 다음과 같이 계산됩니다.

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

그런 다음 정확한 그래디언트 $\left( \nabla_\beta\, \ell(\beta \,;\, \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)}}$를 근사값 $s^{(t)}$로 대체하여 다음을 얻습니다.

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

따라서

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

GLM 사실 도출

이 섹션에서는 이전 섹션에서 사용된 GLM에 대해 자세히 설명하고 결과를 도출합니다. 그런 다음 TensorFlow의 gradients를 사용하여 로그-가능도와 Fisher 정보의 그래디언트에 대한 파생 공식을 수치적으로 검증합니다.

점수 및 Fisher 정보

확률 밀도가 ParseError: KaTeX parse error: Expected '}', got '\right' at position 24: …\cdot | \theta)\̲r̲i̲g̲h̲t̲}_{\theta \in \…인 매개변수 벡터 $\theta$로 매개변수화된 확률 분포군을 고려하세요. 매개변수 벡터 $\theta_0$에서 결과 $y$의 점수는 $y$($\theta_0$에서 평가됨)의 로그-가능도의 그래디언트로 정의됩니다. 즉,

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

주장: 점수의 기대치는 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

주장(Fisher 정보): 점수의 분산은 로그-가능도의 음수 기대 Hessian과 같습니다.

강하지 않은 규칙성 조건(적분에서 미분을 전달할 수 있음)에서는 다음과 같습니다.

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

여기서 $\nabla_\theta^2 F$는 $(i, j)$ 항목이 $\frac{\partial^2 F}{\partial \theta_i \partial \theta_j}$인 Hessian 행렬을 나타냅니다.

이 수식의 좌변은 매개변수 벡터 $\theta_0$에서 $\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T}}$ 분포군의 Fisher 정보라고 합니다.

주장 증명

다음을 가지고 있습니다.

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

$y$에 대한 적분 하에서 $\theta$에 대한 미분을 전달할 수 있도록 하는 강하지 않은 규칙성 조건을 충족하며, 그러면 다음과 같이 됩니다.

$$\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$의 1차 도함수와 2차 도함수입니다.)

증명

이 분포군에 대해 $\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*}$$

과분산 지수군

(스칼라) 과분산 지수군은 밀도가 다음 형식을 취하는 분포군입니다.

$$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$ 함수는 모든 $\theta$에서 $\int p_{\text{OEF}(m, T)}(y\ |\ \theta, \phi=\phi_0), dy = 1$라는 제약 조건에 의해 완전히 결정됩니다. $\phi_0$의 다른 값으로 생성된 $A$는 모두 동일해야 하며, 이에 따라 $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$의 벡터와 연관됩니다. 분포는 과분산 지수군의 구성원이고 매개변수 $\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은 연결 함수와 분포군의 이름을 지정하여 설명합니다(예: "Bernoulli 분포 및 로짓 연결 함수가 있는 GLM"(로지스틱 회귀 모델이라고도 함)). GLM을 완전히 특성화하려면 $h$ 함수도 지정해야 합니다. $h$가 ID이면 $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$를 데이터세트에 피팅하는 데 매우 적합합니다. Fisher 점수와 같은 준뉴턴 방법은 로그 가능도의 그래디언트와 Fisher 정보에 의존하며, 이제 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.$$

Hessian

두 번째 미분하면 곱의 법칙에 의해 다음을 얻습니다.

$$\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를 사용하여 수치적으로 검증하고, Fisher 정보 수식을 tf.hessians를 사용하여 Monte Carlo 추정으로 검증합니다.

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]: Guo-Xun Yuan, Chia-Hua Ho and Chih-Jen Lin. An Improved GLMNET for L1-regularized Logistic Regression. Journal of Machine Learning Research, 13, 2012. http://www.jmlr.org/papers/volume13/yuan12a/yuan12a.pdf

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

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

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