1. 项目概述：为什么我们要手写逻辑回归？

如果你正在学习机器学习，或者想从零开始理解一个分类模型是如何“思考”的，那么“手写一个逻辑回归分类器”几乎是必经之路。这听起来像是一个教科书式的练习，但它的价值远超你的想象。市面上有无数现成的库，比如 scikit-learn 里的 LogisticRegression ，三行代码就能搞定一个分类任务。那我们为什么还要费劲去手写呢？原因很简单： 知其然，更要知其所以然 。当你亲手实现从Sigmoid函数计算、损失函数定义，到用梯度下降或牛顿法迭代更新参数的全过程时，那些原本抽象的概念——比如“概率”、“最大似然”、“优化”——会瞬间变得具体而清晰。这不仅是巩固数学基础，更是培养你调试模型、理解算法收敛性和性能瓶颈的底层能力。今天，我们就抛开框架，从最根本的数学原理出发，一步步构建一个完整的、可用的概率分类器。

2. 核心原理拆解：逻辑回归的“逻辑”是什么？

逻辑回归虽然名字里带“回归”，但它是不折不扣的分类算法，而且是处理二分类问题的利器。它的核心思想非常直观： 不是直接预测类别标签（0或1），而是预测样本属于正类的概率 。这个概率值介于0和1之间，逻辑回归通过一个巧妙的函数将线性回归的无限范围输出映射到这个概率区间，这个函数就是Sigmoid。

2.1 Sigmoid函数：从线性到概率的桥梁

线性回归的假设是 z = w^T * x + b ，其中 w 是权重向量， b 是偏置项， z 的值域是 (-∞, +∞) 。这显然不适合表示概率。Sigmoid函数登场了，其公式为：

σ(z) = 1 / (1 + e^{-z})

这个函数有什么魔力？首先，无论 z 多大或多小， σ(z) 的输出都被压缩在(0, 1)之间，完美符合概率的定义。其次，它是一个单调递增的平滑函数，这意味着 z 越大，属于正类的概率就越高，符合直觉。最后，它的导数有一个非常优美的形式： σ'(z) = σ(z) * (1 - σ(z)) ，这个特性在后续的梯度计算中会大大简化我们的工作。

我们可以这样理解：逻辑回归模型实际上是 P(y=1|x) = σ(w^T * x + b) 。模型通过学习参数 w 和 b ，使得对于正类样本，线性组合 z 尽可能大，从而 σ(z) 接近1；对于负类样本， z 尽可能小， σ(z) 接近0。

2.2 损失函数：交叉熵损失为何是唯一选择？

定义了模型如何输出概率后，我们需要一个标准来衡量模型预测的好坏，这就是损失函数。在逻辑回归中，我们几乎总是使用 二元交叉熵损失 。为什么不用均方误差（MSE）呢？这是新手常有的困惑。

从理论上讲，MSE用于逻辑回归会导致损失函数非凸，存在很多局部极小值，使得优化过程变得困难。而交叉熵损失则是凸函数，能保证我们找到全局最优解（或接近最优的解）。

从直观上理解，交叉熵衡量的是两个概率分布（真实标签分布和预测概率分布）之间的差异。对于单个样本 (x_i, y_i) ，其损失为：

L(y_i, ŷ_i) = -[y_i * log(ŷ_i) + (1 - y_i) * log(1 - ŷ_i)]

其中 ŷ_i = σ(z_i) 是模型预测的概率。这个公式非常巧妙：

• 当真实标签 y_i=1 时，损失变为 -log(ŷ_i) 。如果模型预测概率 ŷ_i 接近1， -log(ŷ_i) 接近0（损失小）；如果预测概率 ŷ_i 接近0， -log(ŷ_i) 会变得非常大（损失大），从而严厉惩罚模型的错误。
• 当真实标签 y_i=0 时，损失变为 -log(1 - ŷ_i) ，逻辑同理。

整个训练集上的损失（成本函数）就是所有样本损失的平均： J(w, b) = (1/m) * Σ L(y_i, ŷ_i) 。我们的优化目标就是找到一组参数 (w, b) ，使得 J(w, b) 最小化。

3. 参数优化：梯度下降与牛顿法的实战抉择

有了损失函数，接下来就是如何找到使其最小化的参数。这是优化算法的战场，我们重点探讨两种最经典的方法：梯度下降和牛顿法。

3.1 梯度下降：稳扎稳打的迭代之道

梯度下降的思想很朴素：沿着当前点损失函数下降最快的方向（负梯度方向）走一小步，不断重复，直到收敛。对于逻辑回归，我们需要计算损失函数 J 关于参数 w 和 b 的梯度。

经过推导（这里涉及对Sigmoid和交叉熵求导，利用链式法则），我们可以得到非常简洁的梯度公式：

• 对于权重 w 的梯度： ∂J/∂w = (1/m) * X^T * (Ŷ - Y)
• 对于偏置 b 的梯度： ∂J/∂b = (1/m) * Σ (ŷ_i - y_i)

其中 X 是特征矩阵， Y 是真实标签向量， Ŷ 是预测概率向量。这个形式是不是很眼熟？它和线性回归的梯度形式在表面上完全一致，但内涵不同，因为这里的 Ŷ 是通过Sigmoid函数计算出来的。

参数更新公式为： w := w - α * ∂J/∂w b := b - α * ∂J/∂b

这里的 α 就是学习率，它是梯度下降中最重要的超参数。学习率太大，可能会在最小值附近震荡甚至发散；学习率太小，收敛速度会慢得令人难以忍受。通常需要从一个较小的值（如0.01）开始尝试，并根据损失曲线进行调整。

实操心得 ：在实现时，一定要将输入特征进行标准化（如Z-score标准化）。这不仅能让梯度下降更快收敛，还能让学习率的选择变得更稳定。想象一下，如果特征A的范围是[0, 1]，特征B的范围是[0, 10000]，那么权重 w 的更新步伐会被特征B主导，导致优化路径扭曲。

3.2 牛顿法：二阶优化的降维打击

梯度下降只利用了一阶导数（梯度）信息，相当于只知道了“下山最陡的方向”。而牛顿法利用了二阶导数（海森矩阵）信息，相当于不仅知道最陡方向，还知道了“地形的曲率”，从而能预测出更优的步长和方向，实现更快的收敛。

对于逻辑回归，牛顿法的参数更新公式为： θ := θ - H^{-1} * ∇J

其中 θ 代表所有参数 [w; b] ， ∇J 是梯度向量， H 是海森矩阵（损失函数 J 关于 θ 的二阶导数矩阵）。牛顿法的核心优势在于它没有学习率这个超参数，其步长由海森矩阵的逆自动决定。在接近最优点时，它通常能实现二次收敛（误差平方级减少），速度远快于梯度下降。

但是，天下没有免费的午餐。牛顿法有两个显著的缺点：

1. 计算成本高 ：海森矩阵的维度是 (n+1) x (n+1) （n是特征数），计算它及其逆矩阵的复杂度是 O(n^3) 。当特征数量很大时（比如上万维），计算将变得不可行。
2. 存储成本高 ：需要存储一个 n x n 的矩阵，内存消耗大。

因此，在实践中：

• 对于特征数不多（例如几百个）的小型数据集，牛顿法通常是首选，因为它收敛迭代次数少，总体时间可能更优。
• 对于高维数据或大数据集，梯度下降（或其变种如随机梯度下降、小批量梯度下降）仍然是主流。

注意事项 ：牛顿法要求海森矩阵必须是正定的，否则更新方向可能不是下降方向。在逻辑回归中，如果数据不是线性可分的，或者特征间存在多重共线性，海森矩阵可能奇异或非正定。在实际实现中，常会加入一个很小的正则项（如 λ * I ）来保证矩阵的可逆性，这实际上等价于使用了L2正则化的逻辑回归。

4. 从零实现：代码细节与避坑指南

理论说得再多，不如一行代码。让我们抛开 sklearn ，用 NumPy 从头构建一个逻辑回归类。这里我将以梯度下降法为例，并指出实现中的关键细节。

4.1 核心类结构设计

首先，我们规划一下类应该有哪些方法：

• __init__ : 初始化参数（权重、偏置）和超参数（学习率、迭代次数）。
• _sigmoid : 静态方法，计算Sigmoid函数。
• _compute_gradient : 根据当前参数计算梯度和损失。
• fit : 训练方法，执行梯度下降迭代。
• predict_proba : 输出预测概率。
• predict : 输出预测类别（默认以0.5为阈值）。

import numpy as np

class LogisticRegressionGD:
    """使用梯度下降法实现逻辑回归。"""
    
    def __init__(self, learning_rate=0.01, n_iters=1000, fit_intercept=True):
        self.lr = learning_rate
        self.n_iters = n_iters
        self.fit_intercept = fit_intercept
        self.weights = None
        self.bias = None
        self.loss_history = []  # 记录损失历史，用于可视化
        
    def _sigmoid(self, z):
        """计算Sigmoid函数，增加数值稳定性处理。"""
        # 防止z过大或过小导致溢出
        z = np.clip(z, -500, 500)
        return 1 / (1 + np.exp(-z))
    
    def _compute_loss(self, y, y_hat):
        """计算交叉熵损失。"""
        # 防止log(0)出现无穷大
        eps = 1e-15
        y_hat = np.clip(y_hat, eps, 1 - eps)
        return -np.mean(y * np.log(y_hat) + (1 - y) * np.log(1 - y_hat))
    
    def fit(self, X, y):
        """
        训练模型。
        参数:
        X: 特征矩阵，形状 (m_samples, n_features)
        y: 标签向量，形状 (m_samples,)
        """
        # 1. 预处理：添加偏置项
        if self.fit_intercept:
            X = np.column_stack([np.ones(X.shape[0]), X])  # 添加一列1
            
        m, n = X.shape
        self.weights = np.zeros(n)  # 初始化参数
        
        # 2. 梯度下降迭代
        for i in range(self.n_iters):
            # 线性组合
            linear_model = np.dot(X, self.weights)
            # 通过Sigmoid得到概率
            y_hat = self._sigmoid(linear_model)
            
            # 计算梯度
            error = y_hat - y
            gradient = np.dot(X.T, error) / m
            
            # 更新参数
            self.weights -= self.lr * gradient
            
            # 记录损失
            loss = self._compute_loss(y, y_hat)
            self.loss_history.append(loss)
            
            # 可选：每100次迭代打印一次损失
            if i % 100 == 0:
                print(f"Iteration {i}: loss = {loss:.4f}")
                
        # 将权重分解回w和b
        if self.fit_intercept:
            self.bias = self.weights[0]
            self.weights = self.weights[1:]
        else:
            self.bias = 0.0
            
        return self
    
    def predict_proba(self, X):
        """预测属于正类的概率。"""
        if self.fit_intercept:
            X = np.column_stack([np.ones(X.shape[0]), X])
            linear_model = np.dot(X, np.concatenate([[self.bias], self.weights]))
        else:
            linear_model = np.dot(X, self.weights)
        return self._sigmoid(linear_model)
    
    def predict(self, X, threshold=0.5):
        """根据阈值将概率转换为类别标签。"""
        proba = self.predict_proba(X)
        return (proba >= threshold).astype(int)

4.2 实现牛顿法版本的关键点

如果你想挑战自己，实现牛顿法版本，核心在于计算海森矩阵 H 。对于逻辑回归，海森矩阵有一个很好的性质： H = (1/m) * X^T * D * X ，其中 D 是一个对角矩阵，其对角线元素 D_ii = ŷ_i * (1 - ŷ_i) 。这是因为每个样本的二阶导数只依赖于其自身的预测值。

class LogisticRegressionNewton:
    """使用牛顿法实现逻辑回归。"""
    
    def __init__(self, n_iters=10, tol=1e-4, fit_intercept=True, reg_lambda=1e-4):
        self.n_iters = n_iters  # 牛顿法迭代次数通常很少
        self.tol = tol  # 收敛容忍度
        self.fit_intercept = fit_intercept
        self.reg_lambda = reg_lambda  # 正则化系数，防止海森矩阵奇异
        self.weights = None
        
    def fit(self, X, y):
        if self.fit_intercept:
            X = np.column_stack([np.ones(X.shape[0]), X])
            
        m, n = X.shape
        self.weights = np.zeros(n)
        
        for i in range(self.n_iters):
            linear_model = np.dot(X, self.weights)
            y_hat = 1 / (1 + np.exp(-linear_model))
            
            # 梯度
            gradient = np.dot(X.T, (y_hat - y)) / m
            
            # 海森矩阵
            D = np.diag(y_hat * (1 - y_hat))
            H = np.dot(X.T, np.dot(D, X)) / m
            # 添加正则项确保可逆
            H += self.reg_lambda * np.eye(n)
            
            # 牛顿更新：求解 H * delta = gradient
            try:
                delta = np.linalg.solve(H, gradient)
            except np.linalg.LinAlgError:
                # 如果求解失败，使用伪逆
                delta = np.dot(np.linalg.pinv(H), gradient)
                
            self.weights -= delta
            
            # 检查收敛：如果参数变化很小则停止
            if np.linalg.norm(delta) < self.tol:
                print(f"Converged at iteration {i}")
                break
                
        if self.fit_intercept:
            self.bias = self.weights[0]
            self.weights = self.weights[1:]
        else:
            self.bias = 0.0
            
        return self
    # ... predict_proba和predict方法同上

踩坑实录 ：在实现牛顿法时，我最初没有添加正则项 reg_lambda ，结果在某个数据集上运行时直接抛出了 LinAlgError: Singular matrix 异常。这是因为当特征存在线性相关，或者某些预测概率接近0或1时，矩阵 D 的对角线元素会非常小，导致海森矩阵 H 接近奇异（不可逆）。加入一个很小的正则项（ λ * I ）是解决这个问题的标准做法，它不仅在数学上稳定了求逆过程，也等价于给模型增加了微弱的L2正则化，防止过拟合。

5. 模型评估与高级话题

模型训练完成后，我们还需要评估其性能，并理解一些进阶概念。

5.1 如何评估你的分类器？

对于二分类问题，准确率（Accuracy）是最直观的指标，但在类别不平衡的数据集上会失灵。更全面的评估需要看混淆矩阵，并计算精确率（Precision）、召回率（Recall）和F1分数。

from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, roc_auc_score

# 假设 y_true 是真实标签， y_pred 是模型预测的类别
accuracy = accuracy_score(y_true, y_pred)
precision = precision_score(y_true, y_pred)
recall = recall_score(y_true, y_pred)
f1 = f1_score(y_true, y_pred)

# 对于概率输出，可以计算AUC-ROC曲线下面积，这是评估概率模型排序能力的黄金标准
y_pred_proba = model.predict_proba(X_test)
roc_auc = roc_auc_score(y_true, y_pred_proba)

绘制ROC曲线能直观展示模型在不同分类阈值下的性能。AUC值越接近1，说明模型区分正负样本的能力越强。

5.2 从二分类到多分类：OvR与Softmax

我们实现的逻辑回归是二分类的。如何扩展到多分类（比如识别手写数字0-9）？有两种主流策略：

1. 一对多（One-vs-Rest, OvR） ：为每个类别训练一个二分类器，将该类视为正类，其余所有类视为负类。预测时，选择输出概率最高的那个分类器对应的类别。这是我们手写逻辑回归最容易扩展的方式。
2. Softmax回归（多项逻辑回归） ：这是逻辑回归在多分类上的直接推广。它将Sigmoid函数替换为Softmax函数，输出一个概率分布，所有类别的概率之和为1。其损失函数也相应变为多类交叉熵损失。Softmax回归的实现比OvR更“原生”，但需要同时优化所有类别的参数。

5.3 正则化：对抗过拟合的武器

当特征很多或样本量相对较少时，模型容易过拟合（在训练集上表现很好，在测试集上表现差）。正则化通过在损失函数中增加一个惩罚项来约束模型参数的大小，从而鼓励模型更简单。

• L1正则化（Lasso） ：在损失函数中加入权重绝对值的和（ λ * ||w||_1 ）。它倾向于产生稀疏解，即让一部分权重直接变为0，从而实现特征选择。
• L2正则化（Ridge） ：在损失函数中加入权重平方和（ λ * ||w||_2^2 ）。它让所有权重都趋近于0，但通常不会等于0，使模型更平滑。

在我们的梯度下降实现中，加入L2正则化非常简单，只需修改梯度计算和损失计算：

# 在损失计算中
loss = self._compute_loss(y, y_hat) + (self.lambda_ / (2*m)) * np.sum(self.weights**2)
# 在梯度计算中
gradient = (np.dot(X.T, error) / m) + (self.lambda_ / m) * self.weights

这里的 lambda_ 是正则化强度超参数，需要交叉验证来确定。

6. 常见问题与调试技巧实录

手写算法时，你会遇到各种预料之外的问题。下面是我在多次实现中总结出的“避坑指南”。

6.1 梯度消失与数值稳定性

问题 ：在计算Sigmoid函数 1/(1+exp(-z)) 时，如果 z 是一个很大的正数， exp(-z) 会下溢为0，导致分母为1，计算结果正确。但如果 z 是一个很大的负数（比如-1000）， exp(-z) 会变成一个天文数字，导致上溢（ overflow ），计算返回 inf 或直接报错。

解决 ：这就是我在 _sigmoid 函数中加入 np.clip(z, -500, 500) 的原因。将 z 的数值范围限制在一个安全区间内，可以彻底避免溢出。这是一种简单粗暴但非常有效的工程化处理。更优雅的做法是分别处理 z 为正和负的情况，但 clip 方法在绝大多数场景下已经足够。

6.2 学习率选择与损失曲线震荡

问题 ：训练时损失不下降，或者像心电图一样剧烈震荡。

诊断与解决 ：

1. 绘制损失曲线 ：这是最重要的调试工具。如果曲线平坦不降，说明学习率可能太小；如果曲线震荡甚至上升，说明学习率太大。
2. 尝试学习率衰减 ：初期使用较大的学习率快速下降，后期使用较小的学习率精细调整。例如： α = α0 / (1 + decay_rate * epoch) 。
3. 使用自适应优化器思想 ：可以手动实现简单的动量（Momentum）来平滑更新方向，减少震荡。更新公式变为： v = β * v + (1-β) * gradient ， w := w - α * v 。其中 β 通常取0.9， v 是速度向量。

6.3 特征尺度不一致导致收敛慢

问题 ：如之前所述，如果特征尺度差异巨大，权重更新会失衡。

解决 ：在调用 fit 方法前，务必对特征进行标准化。最常用的是Z-score标准化： x' = (x - mean(x)) / std(x) 。对于有异常值的数据，也可以使用RobustScaler（使用中位数和四分位距）。这个步骤能显著提升梯度下降的收敛速度和稳定性，对牛顿法也有好处。

6.4 预测概率全部为0.5左右，模型不学习

问题 ：模型输出概率都在0.5附近，没有区分度，准确率接近随机猜测。

诊断 ：

1. 检查数据标签 ：确认你的 y 标签是整数0和1，而不是字符串‘0’和‘1’。
2. 检查梯度计算 ：在代码中打印出前几次迭代的梯度值。如果梯度值非常小（例如 1e-7 ），说明学习信号太弱。可能是特征与标签几乎不相关，或者Sigmoid函数的输入 z 始终在0附近（导致梯度 σ(z)*(1-σ(z)) 很小）。
3. 检查特征工程 ：可能特征本身不具备预测能力，或者需要更复杂的特征组合（非线性变换）。

6.5 牛顿法迭代几次后损失爆炸（NaN）

问题 ：使用牛顿法时，损失突然变成 NaN 。

解决 ：

1. 强制正则化 ：确保海森矩阵 H 添加了正则项 λ * I ， λ 可以设为 1e-4 或 1e-6 。
2. 检查预测概率 ：在计算海森矩阵的对角矩阵 D 时，确保 y_hat 没有非常接近0或1的值（比如 <1e-15 或 >1-1e-15 ），这会导致 D 的对角元素为0，使 H 奇异。在计算 y_hat 后可以加一个裁剪： y_hat = np.clip(y_hat, 1e-15, 1-1e-15) 。
3. 使用更稳定的求解器 ：用 np.linalg.solve 求解线性方程组失败时，可以回退到使用 np.linalg.lstsq （最小二乘解）或 np.linalg.pinv （伪逆），虽然计算慢一些，但更稳健。

手写完一个逻辑回归，你收获的远不止一个可用的分类器。你深入理解了概率映射、最大似然估计、凸优化这些机器学习基石概念的具体运作方式。下次当你轻松调用 model.fit() 时，你会清楚地知道背后发生了什么，参数如何流动，损失如何下降。这种底层的掌控感，是面对更复杂模型和实际生产问题时，进行有效调试和创新的根本。