ヘルプKaggleにTensorFlowグレートバリアリーフを保護チャレンジに参加

一般化線形モデル

TensorFlow.orgで表示GoogleColabで実行GitHubでソースを表示ノートブックをダウンロード

このノートブックでは、実際の例を使用して一般化線形モデルを紹介します。 TensorFlow確率でGLMを効率的にフィッティングするための2つのアルゴリズムを使用して、この例を2つの異なる方法で解決します。密なデータのフィッシャースコアリングと、疎なデータの座標方向の近接勾配降下法です。我々はRの類似の出力に、coordinatewise近位勾配降下の場合に、真係数に近似された係数を比較し、 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\) Aハイパー表す分散液( "分散")、および \(m\)、 \(h\)、 \(T\)、 \(A\) ユーザ指定のモデルによって特徴付けられます家族。

平均 \(Y\) に依存 \(x\) 線形応答の組成によって \(\eta\) と(逆)リンク機能、すなわち:

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

どこ \(g\) 、いわゆるリンク機能です。 TFPでは、リンク機能とモデルの家族の選択は、共同でspecifedさtfp.glm.ExponentialFamilyサブクラス。例は次のとおりです。

TFPは、以上の分布に応じて名前のモデルファミリに好むY以来ではなく、リンク機能tfp.Distribution sがすでにファーストクラスの市民です。場合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]\)、および太字は、これらの機能のベクトルを示しています。これらの期待と分散がどのような分布を超えているかについての完全な詳細は、以下の「GLMファクトの導出」にあります。

このセクションでは、我々は簡単に説明し、TensorFlow確率の2つの組み込みGLMフィッティングアルゴリズム展示:フィッシャースコア( tfp.glm.fit )とcoordinatewise近位勾配降下( 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}\)、左三重$ \({\ textbf {平均値} _T}(\ boldsymbol {\ ETA})、{\ textbf {ヴァー} _T}(\ boldsymbol {\ ETA}を返します)、{\ textbf {Mean} _T} '(\ boldsymbol {\ eta})\ right)$。

私たちはすることをお勧めしますmodelのインスタンスであるtfp.glm.ExponentialFamilyクラス。いくつかの事前に作成された実装が利用可能であるため、ほとんどの一般的なGLMの場合、カスタムコードは必要ありません。

@tf.function(autograph=False)
def fit_model():
  model_coefficients, linear_response, is_converged, num_iter = tfp.glm.fit(
      model_matrix=x, response=y, model=tfp.glm.BernoulliNormalCDF())
  log_likelihood = tfp.glm.BernoulliNormalCDF().log_prob(y, linear_response)
  return (model_coefficients, linear_response, is_converged, num_iter,
          log_likelihood)

[model_coefficients, linear_response, is_converged, num_iter,
 log_likelihood] = [t.numpy() for t in fit_model()]

print(('is_converged: {}\n'
       '    num_iter: {}\n'
       '    accuracy: {}\n'
       '    deviance: {}\n'
       '||w0-w1||_2 / (1+||w0||_2): {}'
      ).format(
    is_converged,
    num_iter,
    np.mean((linear_response > 0.) == y),
    2. * np.mean(log_likelihood),
    np.linalg.norm(model_coefficients_true - model_coefficients, ord=2) /
        (1. + np.linalg.norm(model_coefficients_true, ord=2))
    ))
is_converged: True
    num_iter: 6
    accuracy: 0.75241
    deviance: -0.992436110973
||w0-w1||_2 / (1+||w0||_2): 0.0231555201462

数学的詳細

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

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

対数尤度の勾配の零点を検索するVanillaNewtonの方法は、更新規則に従います。

\[ \beta^{(t+1)}_{\text{Newton} } := \beta^{(t)} - \alpha \left( \nabla^2_\beta\, \ell(\beta\ ;\ \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)} }^{-1} \left( \nabla_\beta\, \ell(\beta\ ;\ \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)} } \]

どこ \(\alpha \in (0, 1]\) ステップサイズを制御するために使用する学習率です。

フィッシャースコアリングでは、ヘッセ行列を負のフィッシャー情報行列に置き換えます。

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

【なお、ここで注意 \(\mathbf{Y} = (Y_i)_{i=1}^{n}\) 、一方、ランダムで \(\mathbf{y}\) 依然として観測された応答のベクトルです。]

以下の「GLMパラメータをデータに適合させる」の式により、これは次のように簡略化されます。

\[ \begin{align*} \beta^{(t+1)} &= \beta^{(t)} + \alpha \left( \mathbf{x}^\top \text{diag}\left( \frac{ \phi\, {\textbf{Mean}_T}'(\mathbf{x} \beta^{(t)})^2 }{ {\textbf{Var}_T}(\mathbf{x} \beta^{(t)}) }\right)\, \mathbf{x} \right)^{-1} \left( \mathbf{x}^\top \text{diag}\left(\frac{ {\textbf{Mean}_T}'(\mathbf{x} \beta^{(t)}) }{ {\textbf{Var}_T}(\mathbf{x} \beta^{(t)}) }\right) \left(\mathbf{T}(\mathbf{y}) - {\textbf{Mean}_T}(\mathbf{x} \beta^{(t)})\right) \right). \end{align*} \]

L1正則化あり

tfp.glm.fit_sparse道具は、GLMはフィッター以上でアルゴリズムに基づいて、まばらなデータセットに適した元、ホーと林2012 。その機能は次のとおりです。

  • L1正則化
  • マトリックス反転なし
  • 勾配とヘッセ行列の評価はほとんどありません。

まず、コードの使用例を示します。アルゴリズムの詳細については、更なる「ためのアルゴリズム詳細に詳述されているtfp.glm.fit_sparse以下」。

model = tfp.glm.Bernoulli()
model_coefficients_start = tf.zeros(x.shape[-1], np.float32)
@tf.function(autograph=False)
def fit_model():
  return tfp.glm.fit_sparse(
    model_matrix=tf.convert_to_tensor(x),
    response=tf.convert_to_tensor(y),
    model=model,
    model_coefficients_start=model_coefficients_start,
    l1_regularizer=800.,
    l2_regularizer=None,
    maximum_iterations=10,
    maximum_full_sweeps_per_iteration=10,
    tolerance=1e-6,
    learning_rate=None)

model_coefficients, is_converged, num_iter = [t.numpy() for t in fit_model()]
coefs_comparison = pd.DataFrame({
  'Learned': model_coefficients,
  'True': model_coefficients_true,
})

print(('is_converged: {}\n'
       '    num_iter: {}\n\n'
       'Coefficients:').format(
    is_converged,
    num_iter))
coefs_comparison
is_converged: True
    num_iter: 1

Coefficients:

学習された係数は、真の係数と同じスパースパターンを持っていることに注意してください。

# Save the learned coefficients to a file.
with tf.io.gfile.GFile('{}/model_coefficients_prox.csv'.format(DATA_DIR), 'w') as f:
  np.savetxt(f, model_coefficients, delimiter=',')

Rさんと比較glmnet

私たちは、Rののそれにcoordinatewise近位勾配降下の出力を比較glmnet同様のアルゴリズムを使用しています、。

注:このセクションを実行するには、Rcolabランタイムに切り替える必要があります。

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\) と \(H\)。

出発点:座標ニュートン法

coordinatewiseニュートン法では、我々はセット \(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)}\)。

ネガティブなフィッシャー情報をヘシアンに置き換えます

さらにバニラ更新ステップのコストを削減するために、我々は設定でき \(H\) かなり正確なヘッセより(以下「データに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*} \]

ここで、 \(r_{\text{L1} } > 0\) 供給される定(L1正則化係数)であり、 \(\text{SoftThreshold}\) によって定義され、ソフトしきい値演算子であります

\[ \text{SoftThreshold}(\beta, \gamma) := \begin{cases} \beta + \gamma &\text{if } \beta < -\gamma \\ 0 &\text{if } -\gamma \leq \beta \leq \gamma \\ \beta - \gamma &\text{if } \beta > \gamma. \end{cases} \]

この更新ルールには、次の2つのインスピレーションを与えるプロパティがあります。これについては以下で説明します。

  1. 限定的な場合に \(r_{\text{L1} } \to 0\) (すなわち、L1正則)、この更新規則は、元の更新規則と同一です。

  2. この更新規則は、不動点がL1正則化最小化問題の解である近接演算子を適用するものとして解釈できます。

$$ \underset{\beta - \beta^{(t)} \in \text{span}{ \text{onehot}(j^{(t)}) } }{\text{arg min} } \left( -\ell(\beta \,;\, \mathbf{x}, \mathbf{y})

  • r_{\text{L1} } \left\lVert \beta \right\rVert_1 \right). $$

縮退ケース \(r_{\text{L1} } = 0\) 回復元の更新ルール

(1)、ノート場合に参照する \(r_{\text{L1} } = 0\) 次いで \(\gamma^{(t)} = 0\)従ってを

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

したがって、

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

不動点が正則化されたMLEである近接演算子

(参照、(2)最初のノートを参照してくださいウィキペディアを任意のためのこと) \(\gamma > 0\)、更新規則

\[ \left(\beta_{\text{exact-prox}, \gamma}^{(t+1)}\right)_{j^{(t)} } := \text{prox}_{\gamma \lVert \cdot \rVert_1} \left( \beta^{(t)}_{j^{(t)} } + \frac{\gamma}{r_{\text{L1} } } \left( \left( \nabla_\beta\, \ell(\beta \,;\, \mathbf{x}, \mathbf{y}) \right)_{\beta = \beta^{(t)} } \right)_{j^{(t)} } \right) \]

満たす(2)、 \(\text{prox}\) 近接演算子である(参照ゆうこの演算子が示され、 \(\mathsf{P}\))。上式の右辺が計算され、ここで

$$

\left(\beta{\text{exact-prox}, \gamma}^{(t+1)}\right){j^{(t)} }

\text{SoftThreshold} \left( \beta^{(t)}{j^{(t)} } + \frac{\gamma}{r{\text{L1} } } \left( \left( \nabla\beta\, \ell(\beta \,;\, \mathbf{x}, \mathbf{y}) \right){\beta = \beta^{(t)} } \right)_{j^{(t)} } ,\ \gamma \right). $$

具体的には、設定\(\gamma = \gamma^{(t)} = -\frac{\alpha\, r_{\text{L1} } }{\left(H^{(t)}\right)_{j^{(t)}, j^{(t)} } }\)(なお、 \(\gamma^{(t)} > 0\) 長い負の対数尤度が凸状である限り)、我々は、更新ルールを得ます

$$

\left(\beta{\text{exact-prox}, \gamma^{(t)} }^{(t+1)}\right){j^{(t)} }

\text{SoftThreshold} \left( \beta^{(t)}{j^{(t)} } - \alpha \frac{ \left( \left( \nabla\beta\, \ell(\beta \,;\, \mathbf{x}, \mathbf{y}) \right){\beta = \beta^{(t)} } \right){j^{(t)} } }{ \left(H^{(t)}\right)_{j^{(t)}, j^{(t)} } } ,\ \gamma^{(t)} \right). $$

我々は、次に、左正確な勾配$ \置き換える(\ナブラ\ベータ\ \エル(\ベータ\; \ \ mathbf {X}、\ mathbf {Y})\右){\ベータ= \ベータ^ {( T)}}その近似と$ \(s^{(t)}\)、取得

{ALIGNを}開始\ \左(\ベータ{\テキスト{正確-PROX} \ガンマ^ {(T)}} ^ {(T + 1)} \右){J ^ {(T)}}&\約\左テキスト{SoftThreshold} \(\ベータ^ {(T)} {J ^ {(T)}} - \アルファ\ FRAC {\(S左^ {(T)} \右){J ^ {( T)}}} {\左(H ^ {(T)} \右){J ^ {(T)}、J ^ {(T)}}} \ \ガンマ^ {(T)} \右) \&= \テキスト{SoftThreshold} \左(\ベータ^ {(T)} {J ^ {(T)}} - \アルファ\、U ^ {(T)}、\ \ガンマ^ {(T)} \正しい)。 \端{ALIGN}

したがって、

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

GLMファクトの導出

このセクションでは、前のセクションで使用されたGLMに関する詳細を説明し、結果を導き出します。その後、我々はTensorFlowの使用gradients数値的に対数尤度とフィッシャー情報のグラデーションの導出式を検証します。

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

パラメータベクトルによってパラメータの確率分布の家族を検討 \(\theta\)確率密度を有する、 \(\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T} }\)。転帰のスコア \(y\) パラメータベクトルで \(\theta_0\) の対数尤度の勾配であると定義される \(y\) (で評価 \(\theta_0\))、すなわち、

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

主張:スコアの期待値はゼロです

穏やかな規則性の条件下で(積分の下で微分を通過させることができます)、

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

証拠

我々は持っています

\[ \begin{align*} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[\text{score}(Y, \theta_0)\right] &:=\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[\left(\nabla_\theta \log p(Y|\theta)\right)_{\theta=\theta_0}\right] \\ &\stackrel{\text{(1)} }{=} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[\frac{\left(\nabla_\theta p(Y|\theta)\right)_{\theta=\theta_0} }{p(Y|\theta=\theta_0)}\right] \\ &\stackrel{\text{(2)} }{=} \int_{\mathcal{Y} } \left[\frac{\left(\nabla_\theta p(y|\theta)\right)_{\theta=\theta_0} }{p(y|\theta=\theta_0)}\right] p(y | \theta=\theta_0)\, dy \\ &= \int_{\mathcal{Y} } \left(\nabla_\theta p(y|\theta)\right)_{\theta=\theta_0}\, dy \\ &\stackrel{\text{(3)} }{=} \left[\nabla_\theta \left(\int_{\mathcal{Y} } p(y|\theta)\, dy\right) \right]_{\theta=\theta_0} \\ &\stackrel{\text{(4)} }{=} \left[\nabla_\theta\, 1 \right]_{\theta=\theta_0} \\ &= 0, \end{align*} \]

ここで、(1)微分の連鎖律、(2)期待値の定義、(3)積分記号の下で微分を渡す(規則性条件を使用)、(4)確率密度の積分は1です。

主張(フィッシャー情報):スコアの分散は、対数尤度の負の期待ヘッセ行列に等しい

穏やかな規則性の条件下で(積分の下で微分を通過させることができます)、

$$ \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \text{score}(Y, \theta_0) \text{score}(Y, \theta_0)^\top

\right]

-\mathbb{E}_{Y \sim p(\cdot | \theta=\theta0)}\left[ \left(\nabla\theta^2 \log p(Y | \theta)\right)_{\theta=\theta_0} \right] $$

ここで、 \(\nabla_\theta^2 F\) そのヘッセ行列表し、 \((i, j)\) エントリである \(\frac{\partial^2 F}{\partial \theta_i \partial \theta_j}\)。

この式の左辺は、家族のフィッシャー情報と呼ばれる \(\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T} }\) パラメータベクトルで \(\theta_0\)。

請求の証明

我々は持っています

\[ \begin{align*} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \left(\nabla_\theta^2 \log p(Y | \theta)\right)_{\theta=\theta_0} \right] &\stackrel{\text{(1)} }{=} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \left(\nabla_\theta^\top \frac{ \nabla_\theta p(Y | \theta) }{ p(Y|\theta) }\right)_{\theta=\theta_0} \right] \\ &\stackrel{\text{(2)} }{=} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \frac{ \left(\nabla^2_\theta p(Y | \theta)\right)_{\theta=\theta_0} }{ p(Y|\theta=\theta_0) } - \left(\frac{ \left(\nabla_\theta\, p(Y|\theta)\right)_{\theta=\theta_0} }{ p(Y|\theta=\theta_0) }\right) \left(\frac{ \left(\nabla_\theta\, p(Y|\theta)\right)_{\theta=\theta_0} }{ p(Y|\theta=\theta_0) }\right)^\top \right] \\ &\stackrel{\text{(3)} }{=} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \frac{ \left(\nabla^2_\theta p(Y | \theta)\right)_{\theta=\theta_0} }{ p(Y|\theta=\theta_0) } - \text{score}(Y, \theta_0) \,\text{score}(Y, \theta_0)^\top \right], \end{align*} \]

ここで、(1)微分のための連鎖律、(2)微分のための商の法則、(3)再び連鎖律を逆に使用しました。

証明を完了するには、それを示すだけで十分です

\[ \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \frac{ \left(\nabla^2_\theta p(Y | \theta)\right)_{\theta=\theta_0} }{ p(Y|\theta=\theta_0) } \right] \stackrel{\text{?} }{=} 0. \]

そのために、積分記号の下で微分を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倍微分は、 \(\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T} }\) によって定義されました

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

満足に対する分化を通過可能に温和規則条件 \(\theta\) に対して一体下 \(y\)次いで、

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

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

(ここで、 \('\) ので、分化を示し \(c'\) と \(c''\) の第一および第二の誘導体である \(c\)。)

証拠

ディストリビューションのこの家族のために、私たちは持っている \(\text{score}(y, \theta_0) = b(y) - c'(\theta_0)\)。最初の式は、その事実から、次の \(\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ \text{score}(y, \theta_0) \right] = 0\)。次に、

\[ \begin{align*} \text{Var}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ b(Y) \right] &= \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ \left(b(Y) - c'(\theta_0)\right)^2 \right] \\ &= \text{the one entry of } \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ \text{score}(y, \theta_0) \text{score}(y, \theta_0)^\top \right] \\ &= \text{the one entry of } -\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ \left(\nabla_\theta^2 \log p(\cdot | \theta)\right)_{\theta=\theta_0} \right] \\ &= -\mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)} \left[ -c''(\theta_0) \right] \\ &= c''(\theta_0). \end{align*} \]

過度に分散した指数型分布族

A(スカラー)は指数ファミリーoverdispersed密度形をとる分布のファミリーであります

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

ここで、 \(m\) と \(T\) 既知のスカラー値関数であり、 \(\theta\) と \(\phi\) スカラーパラメータです。

【ことを注意 \(A\) 過剰決定される:いずれかの \(\phi_0\)、機能 \(A\) 完全に制約によって決定される\(\int p_{\text{OEF}(m, T)}(y\ |\ \theta, \phi=\phi_0)\, dy = 1\)全てに対して \(\theta\)。 \(A\)の異なる値によって生成の \(\phi_0\) すべての機能に制約を置くれ、同じでなければならない \(m\) と \(T\)。]

十分統計量の平均と分散

「対数分配関数の導関数についての補題」と同じ条件で、

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

\right]

A'(\theta) $$

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

\right]

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

証拠

「対数分配関数の導関数についての補題」により、

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

\right]

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

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

\right]

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

結果が期待が(直線的であることから、以下\(\mathbb{E}[aX] = a\mathbb{E}[X]\))及び分散は次数2均質である(\(\text{Var}[aX] = a^2 \,\text{Var}[X]\))。

一般化線形モデル

一般化線形モデルでは、応答変数の予測分布 \(Y\) 観察予測子のベクトルに関連付けられ \(x\)。分布はoverdispersed指数ファミリーのメンバーであり、パラメータ \(\theta\) により置換されて \(h(\eta)\) \(h\) 既知の関数である、 \(\eta := x^\top \beta\) 、いわゆる線形応答であり、 \(\beta\) ベクトルであります学習するパラメータ(回帰係数)の数。一般に、分散パラメータの \(\phi\) あまりにも学ぶことができたが、私たちのセットアップで、私たちは扱います \(\phi\) 知られるように。だから私たちのセットアップは

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

モデル構造は、分布によって特徴付けられる場合 \(p_{\text{OEF}(m, T)}\) 及び機能 \(h\) パラメータに線形応答を変換します。

伝統的に、線形応答からのマッピング \(\eta\) 平均に \(\mu := \mathbb{E}_{Y \sim p_{\text{OEF}(m, T)}(\cdot\, |\, \theta = h(\eta), \phi)}\left[ Y\right]\) 示され、

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

このマッピングは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\) (それぞれ \({\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、Chih-JenLin。 L1正則化ロジスティック回帰のための改善されたGLMNET。機械学習研究のジャーナル、13、2012年http://www.jmlr.org/papers/volume13/yuan12a/yuan12a.pdf

[2]:skd。ソフトしきい値演算子の導出。 2018. https://math.stackexchange.com/q/511106

[3]:ウィキペディアの寄稿者。学習のための近接勾配法。ウィキペディア、フリー百科事典、2018年https://en.wikipedia.org/wiki/Proximal_gradient_methods_for_learning

[4]:Yao-LiangYu。近接演算子。 https://www.cs.cmu.edu/~suvrit/teach/yaoliang_proximity.pdf