廣義線性模型

在 TensorFlow.org 上檢視 在 Google Colab 中執行 在 GitHub 上檢視來源 下載筆記本

在本筆記本中，我們會透過實際範例介紹廣義線性模型。我們會使用兩種不同的方法來解決此範例，亦即使用兩種演算法，以便在 TensorFlow Probability 中有效率地擬合 GLM：適用於密集資料的費雪計分法，以及適用於稀疏資料的座標軸近端梯度下降法。我們會將擬合係數與真實係數進行比較，並在座標軸近端梯度下降法的情況下，與 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 偏好根據 Y 的分配，而非連結函數來命名模型族系，因為 tfp.Distribution 已經是一級公民。如果 tfp.glm.ExponentialFamily 子類別名稱包含第二個字詞，則表示 非典型連結函數。

GLM 有數個顯著的屬性，可有效實作最大概似估計值。這些屬性中最重要的是對數概似 \(\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\) 列是第 \(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 Probability 中的兩種內建 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}). \]

標準牛頓法 (搜尋對數概似梯度的零點) 會遵循更新規則：

\[ \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]\) 是用於控制步幅的學習率。

在費雪計分法中，我們會將 Hessian 替換為負費雪資訊矩陣：

\[ \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 擬合器，其以 Yuan、Ho 和 Lin 2012 中的演算法為基礎。其功能包括：

  • 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 的演算法詳細資訊

我們會將演算法呈現為牛頓法的三種修改。在每一種修改中，\(\beta\) 的更新規則都是以向量 \(s\) 和矩陣 \(H\) 為基礎，其會近似對數概似的梯度和 Hessian。在步驟 \(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\) 的方式不同。

起點：座標軸牛頓法

在座標軸牛頓法中，我們會將 \(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 將梯度近似為一階

\[ \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)}\) 評估的對數概似精確 Hessian。

將負費雪資訊取代為 Hessian

為了進一步降低標準更新步驟的成本，我們可以將 \(H\) 設定為負費雪資訊矩陣 (可使用下方「將 GLM 參數擬合至資料」中的公式有效計算)，而非精確 Hessian：

\[ \begin{align*} H_{\text{Fisher} }^{(t+1)} &:= \mathbb{E}_{Y_i \sim p_{\text{OEF}(m, T)}(\cdot | \theta = h(x_i^\top \beta^{(t+1)}), \phi)} \left[ \left( \nabla^2_\beta\, \ell(\beta\, ;\, \mathbf{x}, \mathbf{Y}) \right)_{\beta = \beta^{(t+1)} } \right] \\ &= -\mathbf{x}^\top \,\text{diag}\left( \frac{ \phi\, {\textbf{Mean}_T}'(\mathbf{x} \beta^{(t+1)})^2 }{ {\textbf{Var}_T}(\mathbf{x} \beta^{(t+1)}) }\right)\, \mathbf{x} \\ s_{\text{Fisher} }^{(t+1)} &:= s_{\text{vanilla} }^{(t+1)} \\ &= \left( \mathbf{x}^\top \,\text{diag}\left(\frac{ {\textbf{Mean}_T}'(\mathbf{x} \beta^{(t+1)}) }{ {\textbf{Var}_T}(\mathbf{x} \beta^{(t+1)}) }\right) \left(\mathbf{T}(\mathbf{y}) - {\textbf{Mean}_T}(\mathbf{x} \beta^{(t+1)})\right) \right) \end{align*} \]

透過近端梯度下降法進行 L1 正規化

如要納入 L1 正規化，我們會將更新規則取代為：

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

更一般的更新規則：

\[ \begin{align*} \gamma^{(t)} &:= -\frac{\alpha\, r_{\text{L1} } }{\left(H^{(t)}\right)_{j^{(t)},\, j^{(t)} } } \\[2mm] \left(\beta_{\text{reg} }^{(t+1)}\right)_j &:= \begin{cases} \beta^{(t+1)}_j &\text{if } j \neq j^{(t)} \\ \text{SoftThreshold} \left( \beta^{(t)}_j - \alpha\, u^{(t)} ,\ \gamma^{(t)} \right) &\text{if } j = j^{(t)} \end{cases} \end{align*} \]

其中 \(r_{\text{L1} } > 0\) 是提供的常數 (L1 正規化係數)，而 \(\text{SoftThreshold}\) 是軟閾值運算子，定義為：

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

此更新規則具有下列兩個啟發性屬性，我們將在下方說明：

  1. 在極限情況 \(r_{\text{L1} } \to 0\) (即無 L1 正規化) 下，此更新規則與原始更新規則相同。

  2. 此更新規則可以解譯為套用鄰近運算子，其固定點是 L1 正規化最小化問題的解：

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

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

退化情況 \(r_{\text{L1} } = 0\) 會復原原始更新規則

如要查看 (1)，請注意，如果 \(r_{\text{L1} } = 0\)，則 \(\gamma^{(t)} = 0\)，因此：

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

因此：

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

鄰近運算子，其固定點是正規化 MLE

如要查看 (2)，首先請注意 (請參閱 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). \)

然後，我們會將精確梯度 $\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，以數值方式驗證對數概似梯度和費雪資訊的推導公式。

分數和費雪資訊

假設一系列機率分配由參數向量 \(\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。

聲明 (費雪資訊)：分數的變異數等於對數概似負期望 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\) 表示 Hessian 矩陣，其 \((i, j)\) 項目為 \(\frac{\partial^2 F}{\partial \theta_i \partial \theta_j}\)。

此方程式的左側稱為族系 \(\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T} }\) 在參數向量 \(\theta_0\) 的**費雪資訊**。

聲明證明

我們有：

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

其中我們使用了 (1) 微分的鏈鎖律、(2) 微分的除法法則、(3) 再次使用鏈鎖律 (反向)。

如要完成證明，只需證明：

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

若要執行此操作，我們會在積分符號下傳遞微分兩次：

\[ \begin{align*} \mathbb{E}_{Y \sim p(\cdot | \theta=\theta_0)}\left[ \frac{ \left(\nabla^2_\theta p(Y | \theta)\right)_{\theta=\theta_0} }{ p(Y|\theta=\theta_0) } \right] &= \int_{\mathcal{Y} } \left[ \frac{ \left(\nabla^2_\theta p(y | \theta)\right)_{\theta=\theta_0} }{ p(y|\theta=\theta_0) } \right] \, p(y | \theta=\theta_0)\, dy \\ &= \int_{\mathcal{Y} } \left(\nabla^2_\theta p(y | \theta)\right)_{\theta=\theta_0} \, dy \\ &= \left[ \nabla_\theta^2 \left( \int_{\mathcal{Y} } p(y | \theta) \, dy \right) \right]_{\theta=\theta_0} \\ &= \left[ \nabla_\theta^2 \, 1 \right]_{\theta=\theta_0} \\ &= 0. \end{align*} \]

關於對數分割函數導數的引理

如果 \(a\)、\(b\) 和 \(c\) 是純量值函數，\(c\) 可微分兩次，使得由下列項目定義的分配族系 \(\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T} }\)：

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

滿足輕微的規律條件 (允許在關於 \(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\) 的一階和二階導數。)

證明

對於此分配族系，我們有 \(\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]\) ) 和變異數的二次齊次性 ( \(\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 的方式是命名其連結函數及其分佈族系，例如「具有白努利分佈和邏輯斯迴歸連結函數的 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{其中}\ \theta = h(x^\top \beta). \end{align*} \]

因此，根據鏈式法則，

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

另外，根據「充分統計量的平均值和變異數」，我們有 \(A'(\theta) = {\text{Mean}_T}(x^\top \beta)\)。因此，根據「聲明：以充分統計量表示 \(h'\)」，我們有

\[ \cdots = \left(T(y) - {\text{Mean}_T}(x^\top \beta)\right) \frac{ {\text{Mean}_T}'(x^\top \beta)}{ {\text{Var}_T}(x^\top \beta)} \,x. \]

黑塞矩陣

再次微分，根據乘積法則，我們得到

\[ \begin{align*} \nabla_\beta^2 \ell(\beta\, ;\, x, y) &= \left[ -A''(h(x^\top \beta))\, h'(x^\top \beta) \right] h'(x^\top \beta)\, x x^\top + \left[ T(y) - A'(h(x^\top \beta)) \right] h''(x^\top \beta)\, xx^\top ] \\ &= \left( -{\text{Mean}_T}'(x^\top \beta)\, h'(x^\top \beta) + \left[T(y) - A'(h(x^\top \beta))\right] \right)\, x x^\top. \end{align*} \]

費雪資訊

根據「充分統計量的平均值和變異數」，我們有

\[ \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 和 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