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 を効率的に適合させるために 2 つのアルゴリズム（フィッシャーのスコア法、疎なデータの場合はスパースな座標近接勾配降下法）を使用して、この例を 2 つの異なる方法で解決します。近似係数を真の係数と比較し、座標近接勾配降下法の場合は、R の同様の glmnet アルゴリズムの出力と比較します。最後に、GLM のいくつかの重要なプロパティの数学的詳細と導出について説明します。

背景情報

一般化線形モデル（GLM）は、変換（リンク関数）にラップされ、指数型分布族に含まれた応答分布がある線形モデル（$\eta = x^\top \beta$）です。リンク機能と応答分布の選択は非常に柔軟なので、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 サブクラス名に 2 番目の単語が含まれている場合、これは非正規リンク関数を示します。

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]$ であり、太字はこれらの関数のベクトル化を示します。これらの期待と分散がどのような分布を超えているかについての完全な詳細は、以下の「一般化線形モデルの詳細の導出」を参照してください。

例

このセクションでは、TensorFlow Probability に組み込まれている 2 つの GLM フィッティングアルゴリズム、フィッシャースコアリング（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 は、フィッシャースコアリングを実装します。これは、引数として以下を取ります。

• 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

数学的詳細

フィッシャースコアリングは、最尤推定値を見つけるための修正されたニュートン法です。

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

対数尤度の勾配 0 を検索する普通のニュートン法は、更新ルールに従います。

$$\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 は、Yuan, Ho and Lin 2012 のアルゴリズムに基づいて、スパースデータセットにより適した GLM フィッターを実装します。機能は次のとおりです。

• 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 と比較する

座標近接勾配降下法の出力を、同様のアルゴリズムを使用する 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 のアルゴリズムの詳細

ニュートン法に対する 3 つの変更のシーケンスとしてアルゴリズムを提示します。それぞれにおいて、$\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*}$$

この更新は、学習率が $\alpha$ のニュートンのようなステップです。最後の部分（L1 正則化）を除いて、以下の変更は、$s$ と $s$ の更新方法のみが異なります。

座標的ニュートン法で始める

座標的ニュートン法では、$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*}$$

勾配とヘッセ行列の評価はほとんどない

対数尤度の勾配とヘッセ行列の計算は高コストになることが多いため、次のように概算することがすすめられます。

• 通常、ヘッセ行列を局所定数として近似し、（近似）ヘッセ行列を使用して勾配を 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)}$ を$\ beta ^ {(t + 1)}$で評価される対数尤度の正確なヘッセ行列に設定します。

負のフィッシャー情報をヘシアンに置き換える

通常の更新ステップのコストをさらに削減するために、正確なヘッセ行列ではなく、負のフィッシャー情報行列（以下の「GLM パラメーターのデータへの適合」の式を使用して効率的に計算できます）に $H$ を設定します。

$$\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{l0}} > 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}$$

この更新ルールには、以下で説明する 2 つのプロパティがあります。

1. 限定的なケース $r_{\text{l0}} \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{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）を確認するには、（Wikipedia を参照）最初に $\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}$ は近接演算子です（Yu を参照してください。この演算子は、$\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).$$

We then replace the exact gradient $\left( \nabla_\beta\, \ell(\beta \,;\, \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)}}$ with its approximation $s^{(t)}$, obtaining

$$\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 に関する詳細を説明し、結果を導き出します。次に、TensorFlow の gradients を使用して、対数尤度とフィッシャー情報の勾配について導出された式を数値的に検証します。

スコアとフィッシャー情報

パラメータベクトル $\theta$ によってパラメータ化され、確率密度 ParseError: KaTeX parse error: Expected '}', got '\right' at position 24: …\cdot | \theta)\̲r̲i̲g̲h̲t̲}_{\theta \in \… を持つ確率分布のファミリをみてみましょう。パラメータベクトル $\theta_0$ での結果 $y$ のスコアは、 $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=\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}$ です。

この方程式の左辺は、パラメータベクトル $\theta_0$ の族 ParseError: KaTeX parse error: Expected '}', got '\right' at position 24: …\cdot | \theta)\̲r̲i̲g̲h̲t̲}_{\theta \in \… のフィッシャー情報量と呼ばれます。

主張の証明

まず、

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

そのために、積分記号の下で微分を 2 回渡します。

$$\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$ は 2 回微分可能であり、分布族 ParseError: KaTeX parse error: Expected '}', got '\right' at position 24: …\cdot | \theta)\̲r̲i̲g̲h̲t̲}_{\theta \in \… は以下によって定義されます。

$$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$ により特徴付けられます。

Traditionally, the mapping from linear response $\eta$ to mean $\mu := \mathbb{E}_{Y \sim p_{\text{OEF}(m, T)}(\cdot\, |\, \theta = h(\eta), \phi)}\left[ Y\right]$ is denoted

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

このマッピングは 1 対 1 である必要があり、その逆関数 $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.$$

ヘッセ行列

積の法則により 2 回微分することにより、以下が得られます。

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

フィッシャー情報

「十分統計量の平均と分散」により、以下が成り立ちます。

$$\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}$ (resp. ${\textbf{Mean}_T}$, resp. ${\textbf{Var}_T}$) がスカラー値関数 $T$ (resp. ${\text{Mean}_T}$, resp. ${\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]: 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