在 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\) 的均值取决于 \(x\)，通过线性响应 \(\eta\) 和（逆）链接函数的组合，即

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

其中 \(g\) 是所谓的链接函数。在 TFP 中，链接函数和模型族的选择由 tfp.glm.ExponentialFamily 子类联合指定。示例包括

• tfp.glm.Normal，也称为“线性回归”
• tfp.glm.Bernoulli，也称为“逻辑回归”
• tfp.glm.Poisson，也称为“泊松回归”
• tfp.glm.BernoulliNormalCDF，也称为“probit 回归”。

TFP 倾向于根据 Y 上的分布来命名模型族，而不是链接函数，因为 tfp.Distribution 已经是第一类公民。如果 tfp.glm.ExponentialFamily 子类名称包含第二个词，则表示 非规范链接函数。

GLM 具有几个显著的特性，这些特性允许对最大似然估计进行有效实现。其中最主要的特性是对对数似然 \(\ell\) 的梯度和费舍尔信息矩阵的简单公式，费舍尔信息矩阵是对数似然负值的 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 拟合算法：Fisher 评分（tfp.glm.fit）和坐标方向近端梯度下降（tfp.glm.fit_sparse）。

合成数据集

让我们假装加载一些训练数据集。

import numpy as np
import pandas as pd
import scipy
import tensorflow.compat.v2 as tf
tf.enable_v2_behavior()
import tensorflow_probability as tfp
tfd = tfp.distributions

def make_dataset(n, d, link, scale=1., dtype=np.float32):
  model_coefficients = tfd.Uniform(
      low=-1., high=np.array(1, dtype)).sample(d, seed=42)
  radius = np.sqrt(2.)
  model_coefficients *= radius / tf.linalg.norm(model_coefficients)
  mask = tf.random.shuffle(tf.range(d)) < int(0.5 * d)
  model_coefficients = tf.where(
      mask, model_coefficients, np.array(0., dtype))
  model_matrix = tfd.Normal(
      loc=0., scale=np.array(1, dtype)).sample([n, d], seed=43)
  scale = tf.convert_to_tensor(scale, dtype)
  linear_response = tf.linalg.matvec(model_matrix, model_coefficients)

  if link == 'linear':
    response = tfd.Normal(loc=linear_response, scale=scale).sample(seed=44)
  elif link == 'probit':
    response = tf.cast(
        tfd.Normal(loc=linear_response, scale=scale).sample(seed=44) > 0,
                   dtype)
  elif link == 'logit':
    response = tfd.Bernoulli(logits=linear_response).sample(seed=44)
  else:
    raise ValueError('unrecognized true link: {}'.format(link))
  return model_matrix, response, model_coefficients, mask

注意：连接到本地运行时。

在本笔记本中，我们使用本地文件在 Python 和 R 内核之间共享数据。为了启用此共享，请使用您有权读取和写入本地文件的同一台机器上的运行时。

x, y, model_coefficients_true, _ = [t.numpy() for t in make_dataset(
    n=int(1e5), d=100, link='probit')]

DATA_DIR = '/tmp/glm_example'
tf.io.gfile.makedirs(DATA_DIR)
with tf.io.gfile.GFile('{}/x.csv'.format(DATA_DIR), 'w') as f:
  np.savetxt(f, x, delimiter=',')
with tf.io.gfile.GFile('{}/y.csv'.format(DATA_DIR), 'w') as f:
  np.savetxt(f, y.astype(np.int32) + 1, delimiter=',', fmt='%d')
with tf.io.gfile.GFile(
    '{}/model_coefficients_true.csv'.format(DATA_DIR), 'w') as f:
  np.savetxt(f, model_coefficients_true, delimiter=',')

无 L1 正则化

函数 tfp.glm.fit 实现 Fisher 评分，它接受一些参数，例如

• model_matrix = \(\mathbf{x}\)
• response = \(\mathbf{y}\)
• model = 可调用函数，它在给定参数 \(\boldsymbol{\eta}\) 时，返回三元组 $\left( {\textbf{Mean}_T}(\boldsymbol{\eta}), {\textbf{Var}_T}(\boldsymbol{\eta}), {\textbf{Mean}_T}'(\boldsymbol{\eta}) \right)$。

我们建议 model 是 tfp.glm.ExponentialFamily 类的实例。有几个预制实现可用，因此对于大多数常见的 GLM，不需要自定义代码。

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

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

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

is_converged: True
    num_iter: 6
    accuracy: 0.75241
    deviance: -0.992436110973
||w0-w1||_2 / (1+||w0||_2): 0.0231555201462

数学细节

Fisher 评分是对牛顿法的一种修改，用于寻找最大似然估计

\[ \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]\) 是一个学习率，用于控制步长。

在 Fisher 评分中，我们用负 Fisher 信息矩阵替换 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)}\) 设置为在 \(\beta^{(t+1)}\) 处评估的对数似然的精确梯度，并将 \(H^{(t+1)}\) 设置为对数似然的精确 Hessian。

用负 Fisher 信息替换 Hessian

为了进一步降低普通更新步骤的成本，我们可以将 \(H\) 设置为负 Fisher 信息矩阵（可以使用下面“将 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_\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} \]

此更新规则具有以下两个启发性的属性，我们将在下面解释

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}\) 是近端算子（参见 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 来数值验证推导的对数似然梯度和 Fisher 信息公式。

得分和 Fisher 信息

考虑一个由参数向量 \(\theta\) 参数化的概率分布族，具有概率密度 \(\left\{p(\cdot | \theta)\right\}_{\theta \in \mathcal{T} }\)。在参数向量 \(\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}\)。

该等式的左侧被称为族 \(\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) \]

满足允许将关于 \(\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*} \]

过度分散指数族

一个（标量）**过度分散指数族** 是一个分布族，其密度函数的形式为

\[ 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\) 确定。由不同 \(\phi_0\) 值产生的 \(A\) 必须全部相同，这给函数 \(m\) 和 \(T\) 施加了一个约束。]

充分统计量的均值和方差

在与“关于对数配分函数导数的引理”相同的条件下，我们有

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

\right]

A'(\theta) $$

以及

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

\right]

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

证明

根据“关于对数配分函数导数的引理”，我们有

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

\right]

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

以及

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

\right]

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

结果然后从期望是线性的 (\(\mathbb{E}[aX] = a\mathbb{E}[X]\)) 以及方差是 2 次齐次的 (\(\text{Var}[aX] = a^2 \,\text{Var}[X]\)) 的事实得出。

广义线性模型

在广义线性模型中，响应变量 \(Y\) 的预测分布与观察到的预测变量向量 \(x\) 相关联。该分布是过度分散指数族的成员，参数 \(\theta\) 被 \(h(\eta)\) 替换，其中 \(h\) 是已知函数，\(\eta := x^\top \beta\) 是所谓的 **线性响应**，\(\beta\) 是要学习的参数向量（回归系数）。一般来说，分散参数 \(\phi\) 也可以学习，但在我们的设置中，我们将 \(\phi\) 视为已知。因此，我们的设置是

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

其中模型结构由分布 \(p_{\text{OEF}(m, T)}\) 和将线性响应转换为参数的函数 \(h\) 来表征。

传统上，从线性响应 \(\eta\) 到均值 \(\mu := \mathbb{E}_{Y \sim p_{\text{OEF}(m, T)}(\cdot\, |\, \theta = h(\eta), \phi)}\left[ Y\right]\) 的映射表示为

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

此映射需要是一对一的，其逆 \(g\) 被称为此 GLM 的 **连接函数**。通常，人们通过命名其连接函数及其分布族来描述 GLM——例如，“具有伯努利分布和 logit 连接函数的 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. \]

黑塞矩阵

第二次微分，根据乘积法则，我们得到

\[ \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}\)（分别为 \({\textbf{Mean}_T}\)，分别为 \({\textbf{Var}_T}\)）表示广播（向量化）函数，该函数将标量值函数 \(T\)（分别为 \({\text{Mean}_T}\)，分别为 \({\text{Var}_T}\)）应用于每个坐标。然后我们有

\[ \begin{align*} \nabla_\beta \ell(\beta\, ;\, \mathbf{x}, \mathbf{y}) &= \sum_{i=1}^{N} \nabla_\beta \ell(\beta\, ;\, x_i, y_i) \\ &= \sum_{i=1}^{N} \left(T(y) - {\text{Mean}_T}(x_i^\top \beta)\right) \frac{ {\text{Mean}_T}'(x_i^\top \beta)}{ {\text{Var}_T}(x_i^\top \beta)} \, x_i \\ &= \mathbf{x}^\top \,\text{diag}\left(\frac{ {\textbf{Mean}_T}'(\mathbf{x} \beta) }{ {\textbf{Var}_T}(\mathbf{x} \beta) }\right) \left(\mathbf{T}(\mathbf{y}) - {\textbf{Mean}_T}(\mathbf{x} \beta)\right) \\ \end{align*} \]

以及

\[ \begin{align*} \mathbb{E}_{Y_i \sim p_{\text{OEF}(m, T)}(\cdot | \theta = h(x_i^\top \beta), \phi)} \left[ \nabla_\beta^2 \ell(\beta\, ;\, \mathbf{x}, \mathbf{Y}) \right] &= \sum_{i=1}^{N} \mathbb{E}_{Y_i \sim p_{\text{OEF}(m, T)}(\cdot | \theta = h(x_i^\top \beta), \phi)} \left[ \nabla_\beta^2 \ell(\beta\, ;\, x_i, Y_i) \right] \\ &= \sum_{i=1}^{N} -\frac{\phi\, {\text{Mean}_T}'(x_i^\top \beta)^2}{ {\text{Var}_T}(x_i^\top \beta)}\, x_i x_i^\top \\ &= -\mathbf{x}^\top \,\text{diag}\left( \frac{ \phi\, {\textbf{Mean}_T}'(\mathbf{x} \beta)^2 }{ {\textbf{Var}_T}(\mathbf{x} \beta) }\right)\, \mathbf{x}, \end{align*} \]

其中分数表示逐元素除法。

数值验证公式

现在，我们使用 tf.gradients 数值验证上述对数似然梯度公式，并使用 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]: 郭训元，何嘉华，林智仁。改进的 GLMNET 用于 L1 正则化逻辑回归。机器学习研究杂志，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]: 俞尧良。近端算子。https://www.cs.cmu.edu/~suvrit/teach/yaoliang_proximity.pdf