线性回归、逻辑回归与优化基础：参考答案#
说明#
以下答案与正文题目逐题对应。封闭推导题给出关键等式；开放编程和实验题给出参考实现、核验标准与结论边界。
给定设计矩阵后，每个观测标签的条件密度为 \((2\pi\sigma^2)^{-1/2}\exp\{-(y_i-\tilde{\bx}_i\trans\cdot\btheta)^2/(2\sigma^2)\}\)。由误差项的独立性，联合似然函数为
\[L(\btheta,\sigma^2) =(2\pi\sigma^2)^{-n/2} \exp\left\{-\frac{1}{2\sigma^2} \sum_{i=1}^n \left(y_i-\tilde{\bx}_i\trans\cdot\btheta\right)^2\right\}.\]因此，对数似然函数为
\[\ell(\btheta,\sigma^2) =-\frac{n}{2}\log(2\pi\sigma^2) -\frac{n}{2\sigma^2}\mathcal{J}(\btheta).\]当 \(\sigma^2>0\) 固定时，第一项与 \(\btheta\) 无关，第二项中 \(n/(2\sigma^2)>0\)，所以最大化 \(\ell(\btheta,\sigma^2)\) 与最小化 \(\mathcal{J}(\btheta)\) 得到相同的 \(\btheta\)。若 \(\sigma^2\) 也未知，则固定 \(\btheta\) 后有 \(\widehat{\sigma}^2(\btheta)=\mathcal{J}(\btheta)\)；当残差平方和为正时，将其代回对数似然所得的剖面对数似然随 \(\mathcal{J}(\btheta)\) 单调递减，因此关于 \(\btheta\) 的最优解仍是最小二乘解。
展开平方损失或直接使用矩阵微分可得
\[\nabla_{\btheta}\mathcal{J}(\btheta) =\frac{2}{n}\bX\trans\cdot (\bX\cdot\btheta-\by), \qquad \boldsymbol{H}_{\mathcal{J}}(\btheta) =\frac{2}{n}\bX\trans\cdot\bX.\]对任意向量 \(\boldsymbol{v}\)，有
\[\boldsymbol{v}\trans\cdot \boldsymbol{H}_{\mathcal{J}}(\btheta)\cdot \boldsymbol{v} =\frac{2}{n} \lVert\bX\cdot\boldsymbol{v}\rVert_2^2 \geq0.\]因此 Hessian 矩阵半正定，\(\mathcal{J}(\btheta)\) 是凸函数。若 \(\bX\) 满列秩，则任意非零 \(\boldsymbol{v}\) 都满足 \(\bX\cdot\boldsymbol{v}\neq\boldsymbol{0}\)，所以 Hessian 矩阵正定，损失函数严格凸并且最小值唯一。若 \(\bX\) 不满列秩，Hessian 矩阵仍然半正定，凸性仍然成立；但零空间中的不同参数可能给出相同预测，因此最小值不一定唯一。
记 \(a_i(\btheta)=\sigma(\tilde{\bx}_i\trans\cdot\btheta)\)。由 \(\sigma'(z)=\sigma(z)\{1-\sigma(z)\}\) 和链式法则可得 \(\partial a_i(\btheta)/\partial\btheta=a_i(\btheta)\{1-a_i(\btheta)\}\tilde{\bx}_i\)。先对平均负对数似然求一次导数，再对梯度中的每一项求导，得到
\[\begin{split}\begin{aligned} \nabla_{\btheta}\mathcal{J}(\btheta) &=\frac{1}{n}\sum_{i=1}^n \{a_i(\btheta)-y_i\}\tilde{\bx}_i,\\ \boldsymbol{H}_{\mathcal{J}}(\btheta) &=\frac{1}{n}\sum_{i=1}^n a_i(\btheta)\{1-a_i(\btheta)\} \tilde{\bx}_i\cdot\tilde{\bx}_i\trans\\ &=\frac{1}{n}\bX\trans\cdot \boldsymbol{R}(\btheta)\cdot\bX, \end{aligned}\end{split}\]其中，\(\boldsymbol{R}(\btheta)\) 是以 \(a_i(\btheta)\{1-a_i(\btheta)\}\) 为第 \(i\) 个对角元素的对角矩阵。对任意向量 \(\boldsymbol{v}\)，有
\[\boldsymbol{v}\trans\cdot \boldsymbol{H}_{\mathcal{J}}(\btheta)\cdot \boldsymbol{v} =\frac{1}{n}\sum_{i=1}^n a_i(\btheta)\{1-a_i(\btheta)\} (\tilde{\bx}_i\trans\cdot\boldsymbol{v})^2 \geq0.\]因此 Hessian 矩阵半正定，负对数似然是凸函数。参数 \(\btheta\) 取有限值时，\(0<a_i(\btheta)<1\)，所以所有权重均为正；若 \(\bX\) 满列秩，则任意非零 \(\boldsymbol{v}\) 至少对应一个非零的 \(\tilde{\bx}_i\trans\cdot\boldsymbol{v}\)，Hessian 矩阵正定，负对数似然严格凸。严格凸只能保证“若有限最小值存在，则它至多有一个”，不能单独保证最小值一定存在。数据完全分离时，可以沿某个参数方向不断增大参数范数，使负对数似然趋近其下确界却不在任何有限参数处达到，因此有限的最大似然估计可能不存在。
交叉熵非负，因为概率在 \((0,1)\) 内时对数不大于0；在标签确定的单样本情形，正确类别概率趋近1时损失下确界为0。 这段论证中的假设不能省略；删除假设时，应用一个最小反例说明结论在哪一步失效。
样本按行时 \(\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}\)。进一步检查：矩阵乘法的内维必须相同，求和轴在结果中消失，广播轴长度只能为1或与目标轴相等，最终梯度必须与对应参数同形。
直接计算 log(sigmoid(z)) 会在线性运算结果的绝对值很大时下溢；把分类阈值固定为0.5也不一定符合代价不对称任务。 最小测试应只改变触发错误的一个条件，并以手算值、维度断言或性质不变量使失败可重复。
可以把截距并入设计矩阵的第一列，只使用
NumPy完成计算。下面的二元交叉熵直接根据线性运算结果计算，避免先求概率后出现 \(\log(0)\)；梯度下降返回每一步的损失，便于检查收敛过程。import numpy as np def sigmoid(z): z = np.asarray(z, dtype=float) out = np.empty_like(z) positive = z >= 0 out[positive] = 1.0 / (1.0 + np.exp(-z[positive])) exp_z = np.exp(z[~positive]) out[~positive] = exp_z / (1.0 + exp_z) return out def binary_cross_entropy_from_logits(y, z): y = np.asarray(y, dtype=float) z = np.asarray(z, dtype=float) if y.shape != z.shape: raise ValueError("y 和 z 的维度必须相同") if not np.all((y == 0) | (y == 1)): raise ValueError("y 只能包含 0 和 1") loss = np.maximum(z, 0.0) + np.log1p(np.exp(-np.abs(z))) - y * z return loss.mean() def logistic_gradient_descent(X, y, learning_rate=0.05, n_steps=5000): X = np.asarray(X, dtype=float) y = np.asarray(y, dtype=float) if X.ndim != 2 or y.ndim != 1 or X.shape[0] != y.size: raise ValueError("X 应为二维数组，且样本数必须与 y 的长度相同") if not np.all(np.isfinite(X)) or not np.all(np.isfinite(y)): raise ValueError("X 和 y 必须只包含有限数值") if not np.all((y == 0) | (y == 1)): raise ValueError("y 只能包含 0 和 1") if learning_rate <= 0 or n_steps <= 0: raise ValueError("学习率和迭代次数必须为正") X_design = np.column_stack((np.ones(X.shape[0]), X)) theta = np.zeros(X_design.shape[1]) loss_history = np.empty(n_steps) for step in range(n_steps): z = X_design @ theta p = sigmoid(z) gradient = X_design.T @ (p - y) / y.size theta -= learning_rate * gradient loss_history[step] = binary_cross_entropy_from_logits( y, X_design @ theta ) return theta, loss_history
核对时，应把“自行实现”和“库函数验证”分开：实现部分只依赖
NumPy；验证部分可以使用scipy.special.expit检查 sigmoid 函数，使用sklearn.metrics.log_loss检查交叉熵，并使用不含正则项的sklearn.linear_model.LogisticRegression检查最终预测概率。示例为from scipy.special import expit from sklearn.metrics import log_loss z_test = np.array([-100.0, -2.0, 0.0, 2.0, 100.0]) y_test = np.array([0.0, 0.0, 1.0, 1.0, 1.0]) np.testing.assert_allclose(sigmoid(z_test), expit(z_test)) np.testing.assert_allclose( binary_cross_entropy_from_logits(y_test, z_test), log_loss(y_test, expit(z_test), labels=[0, 1]), )
比较逻辑回归时，应让两种实现使用相同的数据、截距设置和正则化设置。梯度下降充分收敛后，应重点比较预测概率、交叉熵和分类结果；由于停止条件和求解算法不同，参数不必逐位完全相同。
sigmoid 函数采用正负分支计算，交叉熵根据线性运算结果并使用
logaddexp直接计算；这样避免exp溢出和log(0)。Newton 步用线性方程求解并监控 Hessian 条件数。对logaddexp感兴趣的同学，请参见 logaddexp 的含义与提出背景。记目标成功概率为 \(\pi=\Pr(Y=1)\)。可以选择 \(\pi\in\{0.2,0.5,0.8\}\) 和若干学习率，例如 \(\eta\in\{10^{-3},10^{-2},10^{-1}\}\)。固定特征分布和真实斜率 \(\bw_0\)，然后针对每个 \(\pi\) 调整截距 \(b_0\)，使
\[\frac{1}{n}\sum_{i=1}^{n} \sigma\!\left(b_0+\bx_i\trans\cdot\bw_0\right) \approx\pi,\]再按照 \(Y_i\sim\operatorname{Bernoulli}\{\sigma(b_0+\bx_i\trans\cdot\bw_0)\}\) 生成观测标签。这样只通过截距改变总体成功概率，同时保留特征与观测标签之间的关系。若直接令 \(Y_i\sim\operatorname{Bernoulli}(\pi)\) 且不依赖特征，则真实斜率为零，实验主要检验截距估计，必须在报告中明确说明。
对每个成功概率和随机种子，先生成并保存一份训练集、验证集和测试集；比较不同学习率时必须复用完全相同的数据划分、参数初值和迭代次数。训练过程中保存 \(\btheta^{(t)}\)、训练损失和梯度范数，由此绘制参数更新路径与损失曲线。学习率过小时路径通常移动缓慢，过大时可能振荡或发散；成功概率远离 \(0.5\) 时类别更不平衡，截距也会发生系统变化。
测试集只在训练方案确定后评价一次。除测试交叉熵和准确率外，还应报告 Precision、Recall、F1 或 ROC-AUC，并给出多个固定随机种子下的均值与标准差。准确率会受到成功概率影响，因此不能只凭准确率比较不同 \(\pi\) 下的模型。最终结果表以学习率和成功概率为两个分组变量；其他设置保持不变，才能把参数路径和测试表现的差异归因于这两个研究因素。
每次蒙特卡洛模拟都应先生成一份数据，并让 Newton-Raphson 算法和所有学习率下的梯度下降法复用这份数据、同一个训练集与测试集划分以及同一个参数初值。两种算法分别执行 \(100\) 步参数更新，并在第 \(0,1,\ldots,100\) 步保存参数、训练损失和梯度范数。梯度下降法的更新为
\[\btheta^{(t+1)} =\btheta^{(t)}-eta\nabla\mathcal{J}\!\left(\btheta^{(t)}\right),\]Newton-Raphson 算法则先通过线性方程
\[\boldsymbol{H}_{\mathcal{J}}\!\left(\btheta^{(t)}\right) \cdot\boldsymbol{\Delta}^{(t)} =\nabla\mathcal{J}\!\left(\btheta^{(t)}\right),\]求出 Newton 方向，再令 \(\btheta^{(t+1)}=\btheta^{(t)}-\boldsymbol{\Delta}^{(t)}\)。程序中应使用线性方程求解函数，不要显式计算 Hessian 矩阵的逆。
计时只包括 \(100\) 步参数更新，不包括数据生成、绘图和结果写入。可先预运行一次，再重复计时若干次并报告中位数，同时报告每步平均时间。每次模拟至少记录最终训练损失、测试损失、参数误差 \(\lVert\btheta^{(100)}-\btheta_0\rVert_2\) 和是否出现无穷大或无效数值；多次蒙特卡洛模拟后报告这些量的均值与标准差。还应画出损失曲线和参数更新路径。只比较“都运行了 \(100\) 步”并不完全公平，因此可以同时报告达到同一目标损失所需的步数和时间。
对梯度下降法可选择从小到大的学习率，例如 \(10^{-3}\)、\(10^{-2}\)、\(10^{-1}\) 和 \(1\)。较小学习率通常使损失稳定但下降缓慢；适中的学习率通常能在 \(100\) 步内取得较好的结果；过大的学习率可能使参数在最优点附近振荡，甚至导致损失上升或出现无效数值。这个结论依赖特征尺度，因此比较学习率时必须保持特征预处理方式不变。
Newton-Raphson 算法利用 Hessian 矩阵的曲率信息，在参数维度较低且 Hessian 矩阵条件良好时通常需要较少的更新步数，对学习率也不敏感；但每一步都要构造 Hessian 矩阵并求解线性方程，时间和内存成本较高，在共线性、近完全分离或高维问题中还可能不稳定。梯度下降法每一步只需要梯度，单步计算和存储成本较低，更适合高维或大样本问题；其不足是收敛速度和更新路径对学习率及特征尺度较敏感。最终讨论应以实际记录的时间、损失和参数误差为依据，不能仅根据理论上的迭代次数判断哪种算法更好。