1. 项目缘起：为什么我们要“手写”逻辑回归？

在机器学习入门阶段，逻辑回归（Logistic Regression）几乎是所有人的必经之路。它结构清晰，原理直观，是连接线性模型与概率世界的绝佳桥梁。然而，很多教程和框架（如Scikit-learn）将其封装得过于“黑盒”，一个 LogisticRegression().fit(X, y) 就完成了所有工作。这固然方便，但也让我们错失了深入理解其内部运作机制的机会——比如，损失函数为什么是交叉熵？参数更新除了梯度下降，还有哪些更强大的方法？Sigmoid函数是如何将线性输出转化为概率的？

这次，我们不依赖任何现成的机器学习库，从零开始，用纯Python手写一个完整的逻辑回归分类器。我们的目标不仅仅是“跑通代码”，而是要彻底搞懂从 Sigmoid函数映射 、 极大似然估计推导 ，到 梯度下降（Gradient Descent） 和 牛顿法（Newton‘s Method） 这两种核心优化算法的完整链路。你会发现，亲手实现一遍，那些原本抽象的公式和概念会变得异常具体和牢固。这对于后续理解更复杂的模型（如神经网络）的优化过程，有着不可替代的价值。

2. 逻辑回归的核心：从线性输出到概率预测

逻辑回归的本质，是一个广义线性模型。它的核心思想并不复杂：先用一个线性函数对输入特征进行加权求和，得到一个连续的得分（Logit），然后通过一个非线性函数将这个得分“压缩”到(0, 1)区间内，并将其解释为属于正类的概率。

2.1 Sigmoid函数：概率的“转换器”

这个关键的非线性函数，就是Sigmoid函数，也叫Logistic函数。它的数学形式如下：

[ \sigma(z) = \frac{1}{1 + e^{-z}} ]

其中，( z = w^T x + b )，即我们的线性模型输出。

为什么是Sigmoid？我们可以从几个角度理解：

1. 输出范围完美匹配概率 ：无论输入( z )是多大（正无穷）或多小（负无穷），( \sigma(z) )的输出都被严格限制在(0, 1)之间。这正好符合概率的定义。
2. 良好的数学性质 ：它的导数非常简洁，( \sigma'(z) = \sigma(z)(1 - \sigma(z)) )。这个特性在后续推导梯度时至关重要，能极大简化计算。
3. 可解释性 ：Sigmoid函数是单调递增的。这意味着线性部分( z )越大，模型预测为正类（概率接近1）的信心就越强；( z )越小，预测为负类（概率接近0）的信心就越强。决策边界（通常设为概率0.5）对应着( z = 0 )，即( w^T x + b = 0 )，这是一个超平面。

在实际编码时，我们需要特别注意数值稳定性。当( z )是一个很大的负数时， np.exp(-z) 可能会溢出（变成无穷大）。一个常见的技巧是进行数值裁剪或使用等价公式。例如，我们可以这样实现一个稳定的Sigmoid：

import numpy as np

def sigmoid(z):
    # 数值稳定版本的sigmoid
    # 对于大的正数z， 1/(1+exp(-z)) 直接约等于1，避免exp(-z)下溢为0
    # 对于大的负数z，计算 exp(z) / (1+exp(z))，避免exp(-z)上溢
    return np.where(z >= 0,
                    1 / (1 + np.exp(-z)),
                    np.exp(z) / (1 + np.exp(z)))

2.2 模型假设与决策

有了Sigmoid函数，我们的模型假设就可以写为：

[ P(y=1 | x; w, b) = \sigma(w^T x + b) = \frac{1}{1 + e^{-(w^T x + b)}} ] [ P(y=0 | x; w, b) = 1 - P(y=1 | x; w, b) ]

对于一个新样本( x_{new} )，我们计算( \hat{p} = \sigma(w^T x_{new} + b) )。如果( \hat{p} \ge 0.5 )，我们预测其为正类（y=1），否则为负类（y=0）。这个0.5的阈值可以根据业务需求调整（例如在医疗诊断中，为了降低漏诊率，可能会降低阈值）。

3. 参数学习的基石：极大似然估计与交叉熵损失

模型有了，接下来最关键的问题是：如何找到最优的参数( w )和( b )？我们需要一个准则来衡量参数的好坏，这个准则就是 损失函数（Loss Function） 或 成本函数（Cost Function） 。对于逻辑回归，这个准则源于 极大似然估计（Maximum Likelihood Estimation, MLE） 。

3.1 从极大似然到交叉熵损失

极大似然估计的思想很直观：寻找一组参数，使得在当前参数下，观测到现有这批数据的概率（即似然）最大。

对于单个样本( (x^{(i)}, y^{(i)}) )，其似然可以写为一个紧凑的形式： [ P(y^{(i)} | x^{(i)}; w, b) = (\hat{p}^{(i)})^{y^{(i)}} (1 - \hat{p}^{(i)})^{1 - y^{(i)}} ] 其中，( \hat{p}^{(i)} = \sigma(w^T x^{(i)} + b) )。你可以验证一下，当( y^{(i)}=1 )时，上式等于( \hat{p}^{(i)} )；当( y^{(i)}=0 )时，上式等于( 1 - \hat{p}^{(i)} )。

假设我们有( m )个独立同分布的样本，整个数据集的似然就是所有样本似然的乘积： [ L(w, b) = \prod_{i=1}^{m} (\hat{p}^{(i)})^{y^{(i)}} (1 - \hat{p}^{(i)})^{1 - y^{(i)}} ]

连乘容易导致数值下溢，且不便求导。通常我们取其对数，得到 对数似然（Log-Likelihood） ： [ \ell(w, b) = \log L(w, b) = \sum_{i=1}^{m} \left[ y^{(i)} \log(\hat{p}^{(i)}) + (1 - y^{(i)}) \log(1 - \hat{p}^{(i)}) \right] ]

我们的目标是最大化对数似然( \ell(w, b) )。在优化领域，习惯上我们将最大化问题转化为最小化问题。因此，定义 负对数似然 作为我们的损失函数，再除以样本数( m )得到平均损失，这就是著名的 二元交叉熵损失（Binary Cross-Entropy Loss） ：

[ J(w, b) = -\frac{1}{m} \ell(w, b) = -\frac{1}{m} \sum_{i=1}^{m} \left[ y^{(i)} \log(\hat{p}^{(i)}) + (1 - y^{(i)}) \log(1 - \hat{p}^{(i)}) \right] ]

为什么交叉熵损失是合理的？ 从信息论角度看，交叉熵衡量了真实分布（y）与模型预测分布（p̂）之间的差异。当预测完全正确时（y=1且p̂=1，或y=0且p̂=0），损失为0。当预测完全错误时（y=1但p̂→0），损失会趋近于无穷大，对模型产生巨大的“惩罚”，迫使它更新参数。

3.2 损失函数的代码实现

在实现时，我们同样要关注数值稳定性。 log(0) 是未定义的（ -inf ），而我们的Sigmoid输出理论上不会等于0或1，但数值计算可能非常接近。一个常见的做法是给对数函数输入一个极小的常数（ eps ）进行裁剪。

def binary_cross_entropy_loss(y_true, y_pred):
    """
    计算平均二元交叉熵损失。
    参数:
        y_true: 真实标签，形状 (m, )
        y_pred: 预测概率，形状 (m, )
    返回:
        loss: 标量损失值
    """
    m = y_true.shape[0]
    # 数值稳定：避免log(0)
    eps = 1e-15
    y_pred = np.clip(y_pred, eps, 1 - eps)
    # 计算损失
    loss = -np.mean(y_true * np.log(y_pred) + (1 - y_true) * np.log(1 - y_pred))
    return loss

4. 优化算法一：梯度下降的逐步迭代

有了损失函数( J(w, b) )，我们的任务就变成了一个无约束优化问题：找到使( J )最小的( w )和( b )。梯度下降法是最经典、最直观的一阶优化方法。

4.1 梯度的推导与计算

梯度下降的核心是沿着损失函数梯度的反方向更新参数，因为梯度方向是函数值上升最快的方向，反方向就是下降最快的方向。

我们需要求出损失函数( J )对每个参数( w_j )和( b )的偏导数（即梯度）。推导过程会用到链式法则和Sigmoid函数的导数性质。

令 ( z^{(i)} = w^T x^{(i)} + b )， ( a^{(i)} = \sigma(z^{(i)}) = \hat{p}^{(i)} )。

首先，计算损失对单个样本输出( a^{(i)} )的导数： [ \frac{\partial J}{\partial a^{(i)}} = -\frac{1}{m} \left[ \frac{y^{(i)}}{a^{(i)}} - \frac{1 - y^{(i)}}{1 - a^{(i)}} \right] ]

然后，计算( a^{(i)} )对( z^{(i)} )的导数（Sigmoid导数）： [ \frac{\partial a^{(i)}}{\partial z^{(i)}} = a^{(i)}(1 - a^{(i)}) ]

根据链式法则，损失对( z^{(i)} )的导数为： [ \frac{\partial J}{\partial z^{(i)}} = \frac{\partial J}{\partial a^{(i)}} \cdot \frac{\partial a^{(i)}}{\partial z^{(i)}} = -\frac{1}{m} \left[ \frac{y^{(i)}}{a^{(i)}} - \frac{1 - y^{(i)}}{1 - a^{(i)}} \right] \cdot a^{(i)}(1 - a^{(i)}) ] 化简后，得到一个非常简洁的形式： [ \frac{\partial J}{\partial z^{(i)}} = \frac{1}{m} (a^{(i)} - y^{(i)}) ]

这个结果非常优美！它意味着对于单个样本，损失对线性输出的梯度，就是 预测误差 （预测值 - 真实值）。

接下来，由于( z^{(i)} = w^T x^{(i)} + b = \sum_{j=1}^{n} w_j x_j^{(i)} + b )，我们可以很容易地求出对( w_j )和( b )的梯度： [ \frac{\partial J}{\partial w_j} = \sum_{i=1}^{m} \frac{\partial J}{\partial z^{(i)}} \cdot \frac{\partial z^{(i)}}{\partial w_j} = \frac{1}{m} \sum_{i=1}^{m} (a^{(i)} - y^{(i)}) \cdot x_j^{(i)} ] [ \frac{\partial J}{\partial b} = \sum_{i=1}^{m} \frac{\partial J}{\partial z^{(i)}} \cdot \frac{\partial z^{(i)}}{\partial b} = \frac{1}{m} \sum_{i=1}^{m} (a^{(i)} - y^{(i)}) ]

如果用矩阵形式表示，令 ( X ) 为 ( m \times n ) 的特征矩阵， ( Y ) 和 ( A ) 为 ( m \times 1 ) 的列向量，则： [ \frac{\partial J}{\partial w} = \frac{1}{m} X^T (A - Y) ] [ \frac{\partial J}{\partial b} = \frac{1}{m} \sum (A - Y) \quad \text{(在代码中通常用np.mean处理)} ]

4.2 梯度下降的实现与调参

基于上面的梯度公式，批量梯度下降（Batch Gradient Descent）的更新规则如下： [ w := w - \alpha \cdot \frac{\partial J}{\partial w} ] [ b := b - \alpha \cdot \frac{\partial J}{\partial b} ] 其中，( \alpha ) 是学习率（Learning Rate），是梯度下降中最重要的超参数。

def gradient_descent(X, y, w, b, learning_rate, num_iterations):
    """
    使用批量梯度下降优化逻辑回归参数。
    参数:
        X: 特征矩阵，形状 (m, n)
        y: 标签向量，形状 (m, )
        w: 权重向量，形状 (n, )
        b: 偏置标量
        learning_rate: 学习率
        num_iterations: 迭代次数
    返回:
        w, b, losses: 优化后的参数和损失历史记录
    """
    m = X.shape[0]
    losses = []
    
    for i in range(num_iterations):
        # 1. 前向传播：计算预测概率A
        Z = np.dot(X, w) + b  # 线性部分
        A = sigmoid(Z)        # 概率预测
        
        # 2. 计算损失
        loss = binary_cross_entropy_loss(y, A)
        losses.append(loss)
        
        # 3. 反向传播：计算梯度
        dZ = A - y  # 误差，形状 (m, )
        dw = (1/m) * np.dot(X.T, dZ)  # 梯度 wrt w
        db = (1/m) * np.sum(dZ)       # 梯度 wrt b
        
        # 4. 更新参数
        w = w - learning_rate * dw
        b = b - learning_rate * db
        
        # 可选：每1000次迭代打印一次损失
        if i % 1000 == 0:
            print(f"Iteration {i}: Loss = {loss:.4f}")
    
    return w, b, losses

实操心得与调参技巧：

1. 学习率的选择 ：学习率太大，损失可能会震荡甚至发散；学习率太小，收敛速度会非常慢。一个常用的策略是从一个较大的值（如0.1）开始尝试，观察损失曲线。如果损失震荡，就调小（如0.01， 0.001）；如果下降太慢，可以适当调大。更高级的方法是使用学习率衰减。
2. 特征缩放 ：逻辑回归的损失函数虽然不像线性回归的均方误差那样对特征尺度敏感，但进行特征标准化（Standardization）或归一化（Normalization）通常能加速梯度下降的收敛。因为这样可以让每个特征对梯度的贡献处于同一量级，优化路径更平滑。
3. 迭代停止条件 ：除了固定迭代次数，更合理的做法是设置一个容忍度（tolerance）。当两次迭代之间的损失下降值小于这个容忍度时，就提前停止，认为已经收敛。
4. 初始化 ：参数( w )和( b )通常初始化为0或小的随机数。对于逻辑回归，初始化为0是可行的，因为损失函数是凸的（对于线性可分情况），总能找到全局最优。

5. 优化算法二：牛顿法的二阶收敛

梯度下降是一阶优化方法，它只利用了损失函数的一阶导数（梯度）信息。而牛顿法（Newton‘s Method）是一种二阶优化方法，它同时利用了一阶导数（梯度）和二阶导数（海森矩阵，Hessian Matrix）的信息。在逻辑回归的语境下，牛顿法通常收敛得更快，迭代次数远少于梯度下降。

5.1 牛顿法的原理与优势

牛顿法的核心思想是：在当前位置，用一个二次函数（泰勒二阶展开）来局部近似原损失函数，然后直接找到这个二次函数的最小值点作为下一次的迭代点。

对于我们的损失函数( J(w) )（这里为了简化，将( b )并入( w )，即( w )包含偏置项，( X )增加一列1），其更新公式为： [ w_{new} = w_{old} - H^{-1} \nabla J(w_{old}) ] 其中：

• ( \nabla J(w) ) 是梯度向量（一阶导数）。
• ( H ) 是海森矩阵（Hessian Matrix），即损失函数对( w )的二阶偏导数矩阵。( H_{ij} = \frac{\partial^2 J}{\partial w_i \partial w_j} )。
• ( H^{-1} ) 是海森矩阵的逆。

与梯度下降的更新公式 ( w := w - \alpha \nabla J ) 对比，牛顿法用 ( H^{-1} ) 替代了学习率 ( \alpha )。这个 ( H^{-1} ) 不仅决定了更新的方向，还自动确定了每一步的“最佳”步长。它考虑了损失函数在当前位置的曲率信息：在曲率大的方向（特征值大）迈小步，在曲率小的方向（特征值小）迈大步。

5.2 逻辑回归海森矩阵的推导与计算

对于逻辑回归的交叉熵损失函数，其海森矩阵有一个非常好的性质：它是半正定的（在数据非退化的情况下是正定的），这保证了牛顿法的收敛性。更重要的是，它的形式可以简洁地表示出来。

回顾我们的梯度： [ \nabla J(w) = \frac{1}{m} X^T (A - Y) ]

海森矩阵 ( H ) 是梯度对 ( w ) 的雅可比矩阵。经过推导（过程略），可以得到： [ H = \frac{1}{m} X^T D X ] 其中 ( D ) 是一个 ( m \times m ) 的对角矩阵，其第 ( i ) 个对角线元素为： [ D_{ii} = a^{(i)} (1 - a^{(i)}) ] 这里 ( a^{(i)} = \sigma(w^T x^{(i)}) ) 是第 ( i ) 个样本的预测概率。

这个形式的意义 ：( D_{ii} ) 是Sigmoid函数在 ( z^{(i)} ) 处的导数。当预测概率接近0.5时，( D_{ii} ) 较大（最大为0.25），说明模型在该样本点附近不确定性高，曲率大；当预测概率接近0或1时，( D_{ii} ) 接近0，说明模型很确定，曲率小。海森矩阵通过 ( X^T D X ) 的形式，将每个样本的“不确定性权重”融入了整体的曲率信息中。

5.3 牛顿法的实现与注意事项

牛顿法的实现步骤比梯度下降稍复杂，因为涉及矩阵求逆。

def newtons_method(X, y, w, num_iterations, tol=1e-6):
    """
    使用牛顿法优化逻辑回归参数。
    参数:
        X: 特征矩阵，形状 (m, n)。**注意：这里X应已包含偏置项列（全1列）**。
        y: 标签向量，形状 (m, )
        w: 权重向量（包含偏置），形状 (n, )
        num_iterations: 最大迭代次数
        tol: 收敛容忍度
    返回:
        w, losses: 优化后的参数和损失历史记录
    """
    m = X.shape[0]
    losses = []
    
    for i in range(num_iterations):
        # 1. 前向传播
        Z = np.dot(X, w)
        A = sigmoid(Z)
        
        # 2. 计算损失
        loss = binary_cross_entropy_loss(y, A)
        losses.append(loss)
        
        # 3. 计算梯度
        grad = (1/m) * np.dot(X.T, (A - y))  # 形状 (n, )
        
        # 4. 计算海森矩阵
        # D是对角矩阵，对角线元素为 A * (1 - A)
        D = np.diag(A * (1 - A))  # 形状 (m, m)
        H = (1/m) * np.dot(np.dot(X.T, D), X)  # 形状 (n, n)
        
        # 5. 牛顿法更新: w_new = w_old - H^{-1} * grad
        # 使用np.linalg.solve求解线性方程组 H * delta_w = grad， 比直接求逆更稳定高效
        try:
            delta_w = np.linalg.solve(H, grad)
        except np.linalg.LinAlgError:
            # 如果H是奇异矩阵或接近奇异，添加一个小的正则项（岭回归思想）
            print(f"Iteration {i}: Hessian is singular, adding regularization.")
            H_reg = H + 1e-4 * np.eye(H.shape[0])
            delta_w = np.linalg.solve(H_reg, grad)
        
        w = w - delta_w
        
        # 6. 检查收敛：如果参数变化或梯度范数很小，则停止
        if np.linalg.norm(delta_w) < tol:
            print(f"Converged at iteration {i}.")
            break
            
        if i % 10 == 0: # 牛顿法收敛快，可以少打印几次
            print(f"Iteration {i}: Loss = {loss:.6f}, ||grad|| = {np.linalg.norm(grad):.6f}")
    
    return w, losses

牛顿法的优缺点与实操陷阱：

1. 收敛速度 ：在最优解附近，牛顿法具有二次收敛速率，通常只需不到10次迭代就能达到极高精度，远快于梯度下降的数十上百次迭代。
2. 计算成本 ：每次迭代都需要计算并存储 ( n \times n ) 的海森矩阵及其逆（或解线性方程组），计算复杂度为 ( O(n^3) )。当特征数量 ( n ) 很大（例如上万个）时，这将是不可承受的。因此，牛顿法适用于特征数不太多（几百到几千）的中小规模数据集。
3. 海森矩阵的正定性 ：理论上，逻辑回归的海森矩阵是半正定的。但在实际计算中，由于数值精度或数据共线性，( H ) 可能变成奇异矩阵（不可逆）。代码中我们通过捕获 LinAlgError 并添加一个小的正则项（( H + \lambda I )）来解决，这被称为 正则化牛顿法 。
4. 初始化 ：牛顿法对初始值比梯度下降更敏感。如果初始点离最优点太远，牛顿方向可能不是下降方向。一个稳妥的策略是先用梯度下降跑几轮，得到一个不错的初始点，再切换为牛顿法。

6. 从理论到实践：构建完整的分类器类

现在，我们将前面的所有模块整合起来，构建一个完整的、可复用的逻辑回归分类器类。这个类将支持两种优化算法，并包含预测、评估等完整功能。

import numpy as np

class LogisticRegressionClassifier:
    """
    从零实现的手工逻辑回归分类器。
    """
    def __init__(self, fit_intercept=True, optimizer='gradient_descent', **kwargs):
        """
        初始化分类器。
        参数:
            fit_intercept: 是否拟合偏置项b。
            optimizer: 优化器，可选 'gradient_descent' 或 'newton'。
            **kwargs: 优化器参数，如 learning_rate, num_iterations, tol 等。
        """
        self.fit_intercept = fit_intercept
        self.optimizer = optimizer
        self.params = {} # 存储参数 w, b
        self.loss_history = []
        self.optimizer_kwargs = kwargs
        
    def _add_intercept(self, X):
        """如果需要，为特征矩阵添加一列偏置项（全1）."""
        if self.fit_intercept:
            # 在X的第一列添加一列1
            intercept = np.ones((X.shape[0], 1))
            return np.hstack((intercept, X))
        return X
    
    def _initialize_parameters(self, n_features):
        """初始化模型参数."""
        # 将偏置b作为w向量的第一个元素（如果fit_intercept=True）
        # 这样在牛顿法中处理起来更方便
        if self.fit_intercept:
            # w[0] 将作为偏置b， w[1:] 作为权重
            self.w = np.zeros(n_features + 1)
        else:
            self.w = np.zeros(n_features)
    
    def fit(self, X, y):
        """
        训练逻辑回归模型。
        参数:
            X: 训练特征，形状 (m, n)
            y: 训练标签，形状 (m, )
        """
        m, n = X.shape
        X_processed = self._add_intercept(X)
        self._initialize_parameters(n)
        
        print(f"Training with optimizer: {self.optimizer}")
        if self.optimizer == 'gradient_descent':
            # 从self.w中分离出b和权重
            if self.fit_intercept:
                b_init = self.w[0]
                w_init = self.w[1:]
            else:
                b_init = 0.0
                w_init = self.w
            # 调用梯度下降函数
            learning_rate = self.optimizer_kwargs.get('learning_rate', 0.01)
            num_iterations = self.optimizer_kwargs.get('num_iterations', 10000)
            w_opt, b_opt, losses = gradient_descent(
                X, y, w_init, b_init, learning_rate, num_iterations
            )
            self.loss_history = losses
            if self.fit_intercept:
                self.w = np.concatenate(([b_opt], w_opt))
            else:
                self.w = w_opt
                
        elif self.optimizer == 'newton':
            # 牛顿法：w已包含偏置，X_processed也已包含偏置列
            num_iterations = self.optimizer_kwargs.get('num_iterations', 50)
            tol = self.optimizer_kwargs.get('tol', 1e-6)
            w_opt, losses = newtons_method(
                X_processed, y, self.w, num_iterations, tol
            )
            self.w = w_opt
            self.loss_history = losses
        else:
            raise ValueError(f"Unsupported optimizer: {self.optimizer}")
        
        print("Training completed.")
        return self
    
    def predict_proba(self, X):
        """预测属于正类（y=1）的概率."""
        X_processed = self._add_intercept(X)
        Z = np.dot(X_processed, self.w)
        return sigmoid(Z)
    
    def predict(self, X, threshold=0.5):
        """根据阈值进行类别预测."""
        proba = self.predict_proba(X)
        return (proba >= threshold).astype(int)
    
    def score(self, X, y):
        """计算模型在给定数据上的准确率."""
        y_pred = self.predict(X)
        accuracy = np.mean(y_pred == y)
        return accuracy

使用示例与对比：

# 生成一个简单的二分类数据集
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

# 生成数据
X, y = make_classification(n_samples=1000, n_features=20, n_informative=15, 
                           n_redundant=5, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 特征标准化（对梯度下降很重要）
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# 使用梯度下降训练
print("=== Training with Gradient Descent ===")
lr_gd = LogisticRegressionClassifier(optimizer='gradient_descent', 
                                      learning_rate=0.1, 
                                      num_iterations=5000)
lr_gd.fit(X_train_scaled, y_train)
print(f"Training Accuracy (GD): {lr_gd.score(X_train_scaled, y_train):.4f}")
print(f"Test Accuracy (GD): {lr_gd.score(X_test_scaled, y_test):.4f}")

# 使用牛顿法训练（注意：牛顿法对特征尺度不敏感，但标准化也无害）
print("\n=== Training with Newton's Method ===")
lr_nt = LogisticRegressionClassifier(optimizer='newton', 
                                      num_iterations=20, 
                                      tol=1e-8)
lr_nt.fit(X_train_scaled, y_train)
print(f"Training Accuracy (Newton): {lr_nt.score(X_train_scaled, y_train):.4f}")
print(f"Test Accuracy (Newton): {lr_nt.score(X_test_scaled, y_test):.4f}")

# 对比损失下降曲线
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 5))
plt.plot(lr_gd.loss_history, label='Gradient Descent', alpha=0.7)
plt.plot(lr_nt.loss_history, label="Newton's Method", alpha=0.7)
plt.xlabel('Iteration')
plt.ylabel('Loss (Binary Cross-Entropy)')
plt.title('Loss Convergence Comparison')
plt.legend()
plt.grid(True)
plt.show()

运行这段代码，你会清晰地看到牛顿法在损失收敛速度上的巨大优势，它可能在10次迭代内就达到梯度下降需要数千次迭代才能达到的损失水平。

7. 进阶话题与避坑指南

手写实现让我们能深入控制每一个环节，也必然会遇到一些在调用现成库时被隐藏的问题。

7.1 过拟合与正则化

我们的基本模型没有考虑过拟合。当特征很多或数据有噪声时，模型可能会过于复杂，在训练集上表现很好，但在测试集上表现很差。解决过拟合的标准方法是 正则化（Regularization） 。

在逻辑回归中，最常用的是L2正则化（岭回归），它在损失函数中加入一个惩罚项： [ J_{reg}(w, b) = J(w, b) + \frac{\lambda}{2m} \sum_{j=1}^{n} w_j^2 ] 其中 ( \lambda ) 是正则化强度超参数。注意，通常 不对偏置项 ( b ) 进行正则化 。

加上正则化后，梯度需要相应修改： [ \frac{\partial J_{reg}}{\partial w_j} = \frac{\partial J}{\partial w_j} + \frac{\lambda}{m} w_j ] 矩阵形式为：( \nabla J_{reg} = \nabla J + \frac{\lambda}{m} w )（注意w向量中不包含b）。

海森矩阵也需要修改： [ H_{reg} = H + \frac{\lambda}{m} I ] 其中 ( I ) 是单位矩阵。这实际上使得海森矩阵更加正定，数值稳定性更好。

在代码中如何添加？ 你需要在计算梯度和海森矩阵的步骤中，加上对应的正则化项。 lambda 是一个需要调优的超参数，可以通过交叉验证来选择。

7.2 类别不平衡问题

我们的损失函数是交叉熵，它平等地看待正类和负类样本。但如果数据中正负样本比例悬殊（如99%负样本，1%正样本），模型可能会倾向于将所有样本都预测为多数类，从而得到一个很高的准确率，但这对检测少数类毫无用处。

解决方法：

1. 对损失函数进行加权 ：在交叉熵损失中，给少数类样本更高的权重。 [ J_{weighted} = -\frac{1}{m} \sum_{i=1}^{m} \left[ w_{pos} \cdot y^{(i)} \log(\hat{p}^{(i)}) + w_{neg} \cdot (1 - y^{(i)}) \log(1 - \hat{p}^{(i)}) \right] ] 其中，( w_{pos} ) 和 ( w_{neg} ) 可以根据类别比例设置，例如 ( w_{pos} = \frac{m}{2 \cdot count(y=1)} ), ( w_{neg} = \frac{m}{2 \cdot count(y=0)} )。
2. 重采样 ：对训练数据进行过采样（增加少数类样本）或欠采样（减少多数类样本）。
3. 调整决策阈值 ：不再使用0.5作为阈值。可以根据验证集上的性能（如F1-score、PR曲线）选择一个更优的阈值。

7.3 数值稳定性与计算效率的再思考

1. Log-Sum-Exp Trick ：在计算对数似然或交叉熵时，如果直接计算 log(sigmoid(z)) 或 log(1 - sigmoid(z)) ，在 z 的绝对值很大时可能会遇到数值问题。更稳定的计算方式是： [ \log(\sigma(z)) = -\log(1 + e^{-z}) = z - \log(1 + e^{z}) \quad \text{(当z很大时更稳定)} ] [ \log(1 - \sigma(z)) = -\log(1 + e^{z}) \quad \text{(当z很小时更稳定)} ] 在实际实现损失函数时，可以基于 z 的值选择稳定的计算公式。
2. 牛顿法中的大规模计算 ：对于特征数 n 很大的情况，计算和存储 ( n \times n ) 的海森矩阵及其逆是不现实的。这时可以使用 拟牛顿法（Quasi-Newton Methods） ，如L-BFGS算法。它通过迭代近似海森矩阵的逆，而不需要显式地计算和存储海森矩阵，既能保持较快的收敛速度，又适合大规模问题。Scikit-learn的 LogisticRegression 默认的求解器 lbfgs 就是这一类方法。

7.4 与Scikit-learn的对比验证

作为最终验证，我们可以将自己手写的模型与Scikit-learn官方的实现进行对比，确保我们的实现是正确的。

from sklearn.linear_model import LogisticRegression as SKLogisticRegression

# 使用与我们手写牛顿法类似的配置（L2正则化，无截距，因为我们的X_processed已包含截距）
# 注意：sklearn默认使用L2正则化，C是正则化强度的倒数，C很大相当于正则化很弱。
sk_lr = SKLogisticRegression(penalty='l2', C=1e6, solver='lbfgs', fit_intercept=False, max_iter=1000)
sk_lr.fit(X_train_scaled, y_train)

print("Scikit-learn Model Coefficients (w):", sk_lr.coef_)
print("Our Newton Method Model Coefficients (w):", lr_nt.w[1:]) # 我们的w[0]是偏置
print("Scikit-learn Model Intercept (b):", sk_lr.intercept_)
print("Our Newton Method Model Intercept (b):", lr_nt.w[0])

print(f"\nScikit-learn Test Accuracy: {sk_lr.score(X_test_scaled, y_test):.4f}")
print(f"Our Model Test Accuracy: {lr_nt.score(X_test_scaled, y_test):.4f}")

# 比较预测概率（应该非常接近）
proba_sk = sk_lr.predict_proba(X_test_scaled)[:, 1]
proba_our = lr_nt.predict_proba(X_test_scaled)
print(f"\nMean Absolute Difference in Predicted Probabilities: {np.mean(np.abs(proba_sk - proba_our)):.6f}")

如果一切正确，两个模型的系数应该大致相同，准确率几乎一致，预测概率的差异会非常小（例如小于1e-5）。这个对比是检验我们手写实现正确性的“金标准”。

通过这个从Sigmoid函数、极大似然推导，到梯度下降和牛顿法实现，再到完整分类器构建和进阶话题探讨的完整过程，我们不仅仅是“写出了一个逻辑回归”，更是亲手搭建并理解了概率分类器最核心的优化框架。这份理解，是未来面对更复杂模型时，能够从容拆解其损失函数和优化过程的坚实基础。