从逻辑回归到神经元：参考答案

目录

\[ \begin{align}\begin{aligned}\newcommand{\ba}{\boldsymbol{a}} \newcommand{\bb}{\boldsymbol{b}} \newcommand{\be}{\boldsymbol{e}} \newcommand{\bq}{\boldsymbol{q}} \newcommand{\bk}{\boldsymbol{k}} \newcommand{\bw}{\boldsymbol{w}} \newcommand{\bx}{\boldsymbol{x}} \newcommand{\by}{\boldsymbol{y}} \newcommand{\bz}{\boldsymbol{z}} \newcommand{\bd}{\boldsymbol{d}} \newcommand{\bv}{\boldsymbol{v}} \newcommand{\bs}{\boldsymbol{s}}\\\newcommand{\btheta}{\boldsymbol{\theta}} \newcommand{\bbeta}{\boldsymbol{\beta}} \newcommand{\bgamma}{\boldsymbol{\gamma}} \newcommand{\bsigma}{\boldsymbol{\sigma}} \newcommand{\md}{\mbox{d}} \newcommand{\bmu}{\boldsymbol{\mu}} \newcommand{\bone}{\boldsymbol{1}} \newcommand{\bzero}{\boldsymbol{0}} \newcommand{\bepsilon}{\boldsymbol{\epsilon}} \newcommand{\bphi}{\boldsymbol{\phi}} \newcommand{\bh}{\boldsymbol{h}} \newcommand{\bc}{\boldsymbol{c}} \newcommand{\br}{\boldsymbol{r}} \newcommand{\bQ}{\boldsymbol{Q}} \newcommand{\bK}{\boldsymbol{K}} \newcommand{\bV}{\boldsymbol{V}} \newcommand{\bSigma}{\boldsymbol{\Sigma}} \newcommand{\bg}{\boldsymbol{g}} \newcommand{\bxi}{\boldsymbol{\xi}} \newcommand{\bvarepsilon}{\boldsymbol{\varepsilon}} \newcommand{\bdelta}{\boldsymbol{\delta}} \newcommand{\bq}{\boldsymbol{q}} \newcommand{\bk}{\boldsymbol{k}} \newcommand{\bJ}{\boldsymbol{J}} \newcommand{\bp}{\boldsymbol{p}} \newcommand{\bi}{\boldsymbol{i}} \newcommand{\bo}{\boldsymbol{o}} \newcommand{\bE}{\boldsymbol{E}} \newcommand{\bH}{\boldsymbol{H}} \newcommand{\bL}{\boldsymbol{L}} \newcommand{\bu}{\boldsymbol{u}} \newcommand{\bLambda}{\boldsymbol{\Lambda}} \newcommand{\trans}{^{\rm\scriptsize T}} \newcommand{\var}{\mathrm{var}}\\\newcommand{\bA}{\boldsymbol{A}} \newcommand{\bB}{\boldsymbol{B}} \newcommand{\bC}{\boldsymbol{C}} \newcommand{\bD}{\boldsymbol{D}} \newcommand{\bG}{\boldsymbol{G}} \newcommand{\bI}{\boldsymbol{I}} \newcommand{\bM}{\boldsymbol{M}} \newcommand{\bP}{\boldsymbol{P}} \newcommand{\bS}{\boldsymbol{S}} \newcommand{\bU}{\boldsymbol{U}} \newcommand{\bW}{\boldsymbol{W}} \newcommand{\bX}{\boldsymbol{X}} \newcommand{\bY}{\boldsymbol{Y}} \newcommand{\bZ}{\boldsymbol{Z}} \newcommand{\cotp}{\textcolor[RGB]{48,209,88}{TP}} \newcommand{\cotn}{\textcolor[RGB]{100,210,255}{TN}} \newcommand{\cofp}{\textcolor[RGB]{94,92,230}{FP}} \newcommand{\cofn}{\textcolor[RGB]{191,90,242}{FN}}\\\newcommand{\numcotp}{\textcolor[RGB]{48,209,88}{50}} \newcommand{\numcotn}{\textcolor[RGB]{100,210,255}{30}} \newcommand{\numcofp}{\textcolor[RGB]{94,92,230}{10}} \newcommand{\numcofn}{\textcolor[RGB]{191,90,242}{10}} \DeclareMathOperator*{\argmin}{arg\,min}\end{aligned}\end{align} \]

从逻辑回归到神经元：参考答案#

返回正文练习 · 返回答案索引

说明#

以下答案与正文6道题逐题对应。封闭推导题给出关键等式；开放编程题给出参考实现、核验标准与结论边界。

  1. 神经元先对输入进行仿射变换，再利用激活函数得到输出。仿射变换接收特征和参数，计算可以取任意实数的结果 \(z=\bw\trans\cdot\bx+b\)⁠；sigmoid 函数接收 \(z\)⁠，并将其严格单调地映射到开区间 \((0,1)\)⁠。在逻辑回归中，最终输出 \(a=\sigma(z)\) 可解释为正类条件概率的估计值。这三个概念依次对应计算单元、线性运算和非线性概率映射，不能相互替代。

  2. sigmoid 函数严格单调，并且 \(\sigma(0)=1/2\)⁠，因此

    \[\sigma(z)\geq\frac12 \quad\Longleftrightarrow\quad z\geq 0.\]

    　　当分类规则以 \(0.5\) 为阈值时，两个类别的分界位置满足 \(z=0\)⁠，即 \(\bw\trans\cdot\bx+b=0\)⁠。这个结论依赖于 sigmoid 函数严格单调且阈值恰好为 \(0.5\)⁠；若改用其他阈值，决策边界对应的 \(z\) 也会随之改变。

  3. 对单个样本，\(\bx,\bw\in\mathbb{R}^{d}\)⁠，\(b\)⁠、\(z\) 和 \(a\) 均为标量。对 \(n\) 个按行排列的样本，\(\bX=[\bx_1\trans;\ldots;\bx_n\trans]\in\mathbb{R}^{n\times d}\)⁠、\(\bw\in\mathbb{R}^{d}\)⁠，并且

    \[\bz=\bX\cdot\bw+b\bone\in\mathbb{R}^{n}, \qquad \ba=\sigma(\bz)\in\mathbb{R}^{n}.\]

    　　矩阵 \(\bX\) 的列数必须等于向量 \(\bw\) 的长度；标量 \(b\) 可以广播到所有样本，但程序中应明确检查这一广播只发生在样本轴上。模型参数的梯度必须分别与 \(\bw\) 和 \(b\) 的维度一致。

  4. 单样本实现和批量实现应使用相同的仿射变换与 sigmoid 函数。下面的函数同时返回线性运算结果和概率，便于逐项核验。

    import numpy as np

def sigmoid_stable(z):
    z = np.asarray(z, dtype=float)
    return np.exp(-np.logaddexp(0.0, -z))

def logistic_forward_one(x, w, b):
    x = np.asarray(x, dtype=float)
    w = np.asarray(w, dtype=float)
    if x.ndim != 1 or w.ndim != 1 or x.shape != w.shape:
        raise ValueError("x 和 w 必须是维度相同的一维数组")
    z = x @ w + float(b)
    return z, float(sigmoid_stable(z))

def logistic_forward_batch(X, w, b):
    X = np.asarray(X, dtype=float)
    w = np.asarray(w, dtype=float)
    if X.ndim != 2 or w.ndim != 1 or X.shape[1] != w.size:
        raise ValueError("X 的列数必须等于 w 的长度")
    z = X @ w + float(b)
    return z, sigmoid_stable(z)

X = np.array([[1.0, 2.0], [-1.0, 0.5], [0.0, -2.0]])
w = np.array([0.4, -0.3])
b = 0.2
z_batch, p_batch = logistic_forward_batch(X, w, b)
one_by_one = [logistic_forward_one(x, w, b) for x in X]
np.testing.assert_allclose(z_batch, [v[0] for v in one_by_one])
np.testing.assert_allclose(p_batch, [v[1] for v in one_by_one])

    　　对 logaddexp 感兴趣的同学，请参见 logaddexp 的含义与提出背景⁠。

  5. 为了公平比较循环版和向量化版，应当先生成并保存所有蒙特卡洛样本，再让两个程序复用完全相同的数据、参数初值、学习率和迭代次数。逻辑回归平均损失的梯度为

    \[\nabla\mathcal{J}(\btheta) =\frac1n\widetilde{\bX}\trans\cdot(\bp-\by),\]

    其中，\(\widetilde{\bX}\) 的第一列为1，参数向量 \(\btheta\) 的第一个元素为截距。循环版逐个样本累加同一梯度，向量化版则一次完成矩阵与向量运算。

    from time import perf_counter
import numpy as np

def sigmoid_stable(z):
    z = np.asarray(z, dtype=float)
    return np.exp(-np.logaddexp(0.0, -z))

def generate_samples(repetitions, n, theta_true, seed=1234):
    rng = np.random.default_rng(seed)
    theta_true = np.asarray(theta_true, dtype=float)
    d = theta_true.size - 1
    samples = []
    for _ in range(repetitions):
        X = rng.normal(size=(n, d))
        X_design = np.column_stack((np.ones(n), X))
        p = sigmoid_stable(X_design @ theta_true)
        y = rng.binomial(1, p).astype(float)
        samples.append((X_design, y))
    return samples

def fit_with_loop(X, y, theta_init, learning_rate, steps):
    theta = np.asarray(theta_init, dtype=float).copy()
    n = X.shape[0]
    for _ in range(steps):
        gradient = np.zeros_like(theta)
        for i in range(n):
            p_i = float(sigmoid_stable(X[i] @ theta))
            gradient += (p_i - y[i]) * X[i]
        theta -= learning_rate * gradient / n
    return theta

def fit_vectorized(X, y, theta_init, learning_rate, steps):
    theta = np.asarray(theta_init, dtype=float).copy()
    n = X.shape[0]
    for _ in range(steps):
        p = sigmoid_stable(X @ theta)
        gradient = X.T @ (p - y) / n
        theta -= learning_rate * gradient
    return theta

repetitions = 20
theta_true = np.array([0.5, 1.0, -1.0])
theta_init = np.zeros_like(theta_true)
samples = generate_samples(repetitions, 500, theta_true)

# 预先运行一次，避免把首次调用的额外开销计入正式比较。
fit_with_loop(*samples[0], theta_init, 0.1, 2)
fit_vectorized(*samples[0], theta_init, 0.1, 2)

start = perf_counter()
estimates_loop = np.stack([
    fit_with_loop(X, y, theta_init, 0.1, 100)
    for X, y in samples
])
time_loop = perf_counter() - start

start = perf_counter()
estimates_vectorized = np.stack([
    fit_vectorized(X, y, theta_init, 0.1, 100)
    for X, y in samples
])
time_vectorized = perf_counter() - start

np.testing.assert_allclose(
    estimates_loop, estimates_vectorized, rtol=1e-10, atol=1e-10
)
print("循环版参数均值：", estimates_loop.mean(axis=0))
print("向量化版参数均值：", estimates_vectorized.mean(axis=0))
print("最大参数差异：", np.max(np.abs(
    estimates_loop - estimates_vectorized
)))
print("循环版时间：", time_loop)
print("向量化版时间：", time_vectorized)
print("循环版与向量化版时间之比：", time_loop / time_vectorized)

    　　两个程序实现的是同一个梯度公式，因此在相同数据和设置下，参数估计结果应当在浮点误差范围内一致。向量化版把主要运算交给底层数值计算程序库，通常明显快于在 Python 中逐样本循环；但具体时间比取决于样本量、迭代次数、硬件和数值计算程序库。正式报告应重复计时并报告典型运行时间及波动，不能只依据一次计时下结论。蒙特卡洛参数均值接近真实参数只能说明当前样本量和优化设置下结果合理，不能代替对收敛性与有限样本误差的分析。

  6. sigmoid 函数若直接计算为 \(\{1+\exp(-z)\}^{-1}\)⁠，会在 \(z\) 的绝对值很大时使指数溢出或下溢。可以按 \(z\) 的正负分别计算，也可以使用 logaddexp 得到稳定实现。二元交叉熵应直接根据线性运算结果写成 \(\operatorname{logaddexp}(0,z)-yz\)⁠，避免先计算接近0或1的概率后再取对数。测试时应加入很大的正数、很小的负数和接近0的输入，并确认概率、损失及梯度均为有限值。对 logaddexp 感兴趣的同学，请参见 logaddexp 的含义与提出背景⁠。

目录