从零到一：用Python与NumPy亲手构建线性回归与逻辑回归模型

最近几年，机器学习的热度居高不下，但很多学习者在看完了吴恩达教授等经典课程的理论讲解后，面对空白的代码编辑器，依然会感到无从下手。理论公式和实际可运行的代码之间，似乎总隔着一道鸿沟。这篇文章，就是为你——一位有一定Python基础，渴望将机器学习理论付诸实践的学习者——准备的实战指南。我们将抛开现成的scikit-learn库，仅使用最基础的NumPy，从零开始，一步步推导并实现线性回归和逻辑回归的核心算法。这不仅仅是“复现”，更是一次深入理解模型内部运作机制、掌握向量化编程思维、并学会调试和优化模型的绝佳旅程。你会发现，亲手搭建这些基础模型，对理解更复杂的深度学习框架有着不可替代的价值。

1. 环境搭建与数据准备：为实战打下坚实基础

在开始编写任何机器学习代码之前，一个清晰、可复现的开发环境是必不可少的。我强烈建议使用conda或venv来创建独立的Python虚拟环境，这能有效避免不同项目间的包版本冲突。对于本教程，我们需要的核心库非常简单：NumPy用于数值计算，Matplotlib用于可视化，或许再加上pandas来方便地加载数据。

# 使用 conda 创建环境并安装依赖（推荐）
conda create -n ml_from_scratch python=3.9
conda activate ml_from_scratch
conda install numpy matplotlib pandas jupyter

# 或者使用 pip 和 venv
python -m venv ml_env
source ml_env/bin/activate  # Windows 下使用 ml_env\Scripts\activate
pip install numpy matplotlib pandas

数据是模型的燃料。为了聚焦于算法本身，我们可以自己生成一些模拟数据。例如，对于线性回归，我们可以生成符合 y = 2x + 1 + 噪声 规律的数据点。这样做的好处是，我们预先知道了“标准答案”，可以清晰地验证模型学习的效果。

import numpy as np
import matplotlib.pyplot as plt

# 设置随机种子以保证结果可复现
np.random.seed(42)

# 生成线性回归模拟数据
def generate_linear_data(num_samples=100):
    X = 2 * np.random.rand(num_samples, 1)  # 特征 X，范围在[0, 2)
    # 真实关系: y = 2x + 1 + 高斯噪声
    y = 2 * X + 1 + np.random.randn(num_samples, 1) * 0.5
    return X, y

X_train, y_train = generate_linear_data()
plt.scatter(X_train, y_train)
plt.xlabel('Feature X')
plt.ylabel('Target y')
plt.title('Generated Linear Regression Data')
plt.show()

提示：在机器学习项目中，养成设置随机种子（如 np.random.seed(42)）的习惯至关重要。它能确保每次运行代码时，随机生成的数据、初始化的参数都是一致的，这对于调试和比较不同算法效果非常有帮助。

对于逻辑回归，我们需要生成可用于二分类的数据。scikit-learn的 make_classification 或 make_moons 函数是很好的选择，但为了纯粹性，我们也可以手动构造一个简单的二维数据集，确保两类数据线性可分或具有清晰的决策边界。

2. 线性回归的向量化实现与梯度下降优化

线性回归是理解机器学习优化过程的绝佳起点。其核心思想是找到一组参数（权重 w 和偏置 b），使得模型预测值 ŷ = wX + b 与真实值 y 之间的差距（用损失函数衡量）最小。

2.1 定义模型与损失函数

首先，我们定义假设函数（预测函数）。这里的关键是使用向量化操作，一次性对整个训练集进行计算，这比使用循环效率高出几个数量级。

def linear_hypothesis(X, w, b):
    """
    线性回归假设函数 (向量化实现)
    参数:
        X: 特征矩阵，形状 (m, n)，m为样本数，n为特征数
        w: 权重向量，形状 (n, 1)
        b: 偏置标量
    返回:
        y_pred: 预测值，形状 (m, 1)
    """
    # 核心操作: y_pred = X * w + b
    # 这里使用点积 np.dot，对于单特征，w是标量；多特征时，w是向量
    y_pred = np.dot(X, w) + b
    return y_pred

接下来，我们需要一个标尺来衡量预测的好坏，这就是损失函数（或成本函数）。对于线性回归，最常用的是均方误差（MSE）。

def compute_mse_cost(y_true, y_pred):
    """
    计算均方误差损失 (Mean Squared Error)
    参数:
        y_true: 真实标签，形状 (m, 1)
        y_pred: 模型预测值，形状 (m, 1)
    返回:
        cost: 标量损失值
    """
    m = y_true.shape[0]
    # 向量化计算每个样本的误差平方，然后求和取平均
    cost = (1 / (2 * m)) * np.sum((y_pred - y_true) ** 2)
    return cost

注意：损失函数中的 1/(2*m) 而不是 1/m 是学术上的一个常见写法，其目的是在后续求梯度时，平方项求导产生的系数2能与分母的2约掉，使得梯度表达式更简洁。这对最终优化结果没有影响。

2.2 实现梯度下降算法

梯度下降是优化损失函数的引擎。其原理是计算损失函数关于每个参数的梯度（导数），然后沿着梯度反方向（即下降最快的方向）更新参数。向量化实现能同时更新所有参数。

def batch_gradient_descent(X, y, w, b, learning_rate, num_iterations):
    """
    批量梯度下降算法
    参数:
        X: 特征矩阵 (m, n)
        y: 真实标签 (m, 1)
        w: 初始权重 (n, 1)
        b: 初始偏置 (标量)
        learning_rate: 学习率 α
        num_iterations: 迭代次数
    返回:
        w, b: 优化后的参数
        cost_history: 每次迭代的损失记录，用于绘图
    """
    m = X.shape[0]
    cost_history = []

    for i in range(num_iterations):
        # 1. 计算当前参数的预测值
        y_pred = linear_hypothesis(X, w, b)

        # 2. 计算梯度 (向量化形式)
        # 对w的梯度: dj_dw = (1/m) * X^T * (y_pred - y)
        # 对b的梯度: dj_db = (1/m) * sum(y_pred - y)
        error = y_pred - y
        dj_dw = (1 / m) * np.dot(X.T, error)
        dj_db = (1 / m) * np.sum(error)

        # 3. 同时更新参数
        w = w - learning_rate * dj_dw
        b = b - learning_rate * dj_db

        # 4. 记录当前损失
        cost = compute_mse_cost(y, y_pred)
        cost_history.append(cost)

        # 可选：每1000次迭代打印一次损失
        if i % 1000 == 0:
            print(f"Iteration {i:4d}: Cost {cost:.6f}")

    return w, b, cost_history

现在，让我们在一个简单的单变量数据集上运行这个算法，并观察其学习过程。

# 准备数据 (添加偏置列X0=1，便于将b并入w进行向量化，但这里我们分开处理以清晰理解)
X = X_train  # 形状 (100, 1)
y = y_train

# 初始化参数
w_initial = np.random.randn(1, 1)  # 随机初始化权重
b_initial = 0.0
learning_rate = 0.01
iterations = 5000

# 运行梯度下降
w_final, b_final, costs = batch_gradient_descent(X, y, w_initial, b_initial, learning_rate, iterations)

print(f"\n最终参数: w = {w_final[0,0]:.4f}, b = {b_final:.4f}")
print(f"真实关系近似为: y = {w_final[0,0]:.4f} * x + {b_final:.4f}")

# 绘制学习曲线
plt.plot(range(iterations), costs)
plt.xlabel('Iteration')
plt.ylabel('Cost (MSE)')
plt.title('Gradient Descent: Cost vs. Iterations')
plt.grid(True)
plt.show()

通过绘制损失随迭代次数下降的曲线，我们可以直观判断学习率是否合适、算法是否收敛。一条平滑下降并最终趋于平稳的曲线是健康的标志。

2.3 特征缩放与学习率调优实战

当特征尺度差异巨大时（例如，房屋面积 x1 范围是 [500, 5000]，而房间数量 x2 范围是 [1, 5]），梯度下降会变得非常缓慢，路径曲折。特征缩放是解决这一问题的关键技术，它能将不同特征的值映射到相近的范围内，最常见的方法是Z-score标准化。

缩放方法	公式	优点	缺点
Z-score标准化	x' = (x - μ) / σ	处理后数据均值为0，标准差为1，适用于大多数情况。	需要计算均值和标准差。
Min-Max归一化	x' = (x - min) / (max - min)	将值缩放到 [0, 1] 区间，直观。	对异常值敏感。

def z_score_normalize(X):
    """
    对特征矩阵X进行Z-score标准化 (按列/特征进行)
    参数:
        X: 原始特征矩阵 (m, n)
    返回:
        X_norm: 标准化后的矩阵
        mu: 每个特征的均值 (用于后续对新数据标准化)
        sigma: 每个特征的标准差
    """
    mu = np.mean(X, axis=0)   # 计算每个特征的均值，形状 (1, n)
    sigma = np.std(X, axis=0) # 计算每个特征的标准差
    # 避免除零，给标准差为0的特征一个很小的值
    sigma = np.where(sigma == 0, 1, sigma)
    X_norm = (X - mu) / sigma
    return X_norm, mu, sigma

# 假设我们有一个多特征数据集 X_multi
# X_norm, mu, sigma = z_score_normalize(X_multi)
# 在训练模型时使用 X_norm
# 当有新数据需要预测时，必须使用训练时计算的 mu 和 sigma 进行同样的变换

学习率（α） 是梯度下降中最重要的超参数之一。它控制着每次参数更新的步长。

• α太小：收敛速度极慢，需要大量迭代。
• α太大：可能越过最低点，导致损失震荡甚至发散。

一个实用的调优策略是进行学习率扫描：尝试一系列呈指数增长或衰减的学习率（例如 [0.001, 0.003, 0.01, 0.03, 0.1]），运行少量迭代（如1000次），绘制每个学习率对应的损失下降曲线。选择那个使损失快速且平稳下降的学习率。

def test_learning_rates(X, y, learning_rates):
    """
    测试不同学习率对梯度下降收敛的影响
    """
    plt.figure(figsize=(10, 6))
    for lr in learning_rates:
        w = np.random.randn(X.shape[1], 1)
        b = 0.0
        _, _, costs = batch_gradient_descent(X, y, w, b, lr, num_iterations=1000)
        plt.plot(costs, label=f'lr={lr}')
    plt.xlabel('Iterations')
    plt.ylabel('Cost')
    plt.title('Cost vs. Iterations for Different Learning Rates')
    plt.legend()
    plt.grid(True)
    plt.yscale('log')  # 使用对数坐标更易观察
    plt.show()

# 调用测试
# test_learning_rates(X_norm, y, [0.001, 0.01, 0.1, 0.5])

3. 逻辑回归：从概率视角理解分类

线性回归预测连续值，而逻辑回归解决的是二分类问题（是/否，0/1）。其核心创新在于引入了 Sigmoid函数（或称逻辑函数），将线性组合 z = wX + b 的输出映射到 (0, 1) 区间，并将其解释为样本属于正类的概率。

3.1 Sigmoid函数与决策边界

Sigmoid函数的数学表达式为：g(z) = 1 / (1 + e^{-z})。当 z 很大时，g(z) 接近1；当 z 很小时，g(z) 接近0；在 z=0 时，g(z)=0.5。

def sigmoid(z):
    """
    Sigmoid激活函数
    参数:
        z: 标量或NumPy数组
    返回:
        经过sigmoid变换的值，范围在(0,1)
    """
    # 使用 np.exp 进行向量化计算
    return 1 / (1 + np.exp(-z))

# 可视化Sigmoid函数
z = np.linspace(-10, 10, 100)
g = sigmoid(z)
plt.plot(z, g)
plt.axvline(x=0, color='k', linestyle='--', alpha=0.3)
plt.axhline(y=0.5, color='k', linestyle='--', alpha=0.3)
plt.xlabel('z')
plt.ylabel('g(z)')
plt.title('Sigmoid Function')
plt.grid(True)
plt.show()

逻辑回归的假设函数因此定义为：h(x) = g(wX + b) = sigmoid(wX + b)。我们通常设定一个阈值（默认为0.5）来做出分类决策：

• 如果 h(x) >= 0.5，预测为类别1。
• 如果 h(x) < 0.5，预测为类别0。

这个阈值0.5对应的 z = wX + b = 0，在特征空间中就构成了一条直线（对于二维特征）或一个超平面，称为决策边界。模型学习的目标，就是找到参数 w 和 b，使得这个决策边界能最好地分开两类数据。

3.2 交叉熵损失函数及其梯度

线性回归使用的均方误差损失在逻辑回归中不再适用，因为它会导致损失函数非凸，存在多个局部极小值。逻辑回归使用交叉熵损失函数，它源于最大似然估计的思想，能很好地衡量预测概率分布与真实标签分布的差异。

对于单个样本 (x^(i), y^(i))，其损失为： L(f(x^(i)), y^(i)) = -[y^(i) * log(f(x^(i))) + (1 - y^(i)) * log(1 - f(x^(i)))]

整个训练集的成本函数是所有样本损失的平均值。这个函数是凸的，保证了梯度下降能找到全局最优解。

def compute_logistic_cost(y_true, y_pred_prob):
    """
    计算逻辑回归的交叉熵损失 (向量化实现)
    参数:
        y_true: 真实标签 (0或1)，形状 (m, 1)
        y_pred_prob: 预测为正类的概率，形状 (m, 1)
    返回:
        cost: 标量损失值
    """
    m = y_true.shape[0]
    # 避免log(0)导致数值问题，对预测概率进行裁剪
    epsilon = 1e-15
    y_pred_prob = np.clip(y_pred_prob, epsilon, 1 - epsilon)

    # 交叉熵损失计算
    cost = - (1 / m) * np.sum(y_true * np.log(y_pred_prob) + (1 - y_true) * np.log(1 - y_pred_prob))
    return cost

def logistic_gradient_descent(X, y, w, b, learning_rate, num_iterations):
    """
    逻辑回归的梯度下降算法
    """
    m = X.shape[0]
    cost_history = []

    for i in range(num_iterations):
        # 1. 计算线性部分和概率预测
        z = np.dot(X, w) + b
        y_pred_prob = sigmoid(z)  # f(x)

        # 2. 计算梯度
        error = y_pred_prob - y  # 注意这个形式与线性回归惊人地相似
        dj_dw = (1 / m) * np.dot(X.T, error)
        dj_db = (1 / m) * np.sum(error)

        # 3. 更新参数
        w = w - learning_rate * dj_dw
        b = b - learning_rate * dj_db

        # 4. 记录损失
        cost = compute_logistic_cost(y, y_pred_prob)
        cost_history.append(cost)

        if i % 1000 == 0:
            print(f"Iteration {i:4d}: Cost {cost:.6f}")

    return w, b, cost_history

注意：在计算交叉熵损失时，我们使用 np.clip 将预测概率限制在 [epsilon, 1-epsilon] 范围内。这是因为当预测概率非常接近0或1时，log(0) 会导致数值计算问题（负无穷）。这是一个非常重要的实战技巧。

3.3 模型评估与决策边界可视化

训练完逻辑回归模型后，我们需要评估其性能。对于二分类，最基本的评估指标是准确率。

def predict(X, w, b, threshold=0.5):
    """
    使用训练好的逻辑回归模型进行预测
    参数:
        X: 特征矩阵
        w, b: 模型参数
        threshold: 分类阈值
    返回:
        y_pred_class: 预测的类别 (0或1)
        y_pred_prob: 预测为正类的概率
    """
    z = np.dot(X, w) + b
    y_pred_prob = sigmoid(z)
    y_pred_class = (y_pred_prob >= threshold).astype(int)
    return y_pred_class, y_pred_prob

def accuracy_score(y_true, y_pred):
    """计算分类准确率"""
    return np.mean(y_true == y_pred)

# 假设 X_test, y_test 是测试集
# y_pred_class, _ = predict(X_test, w_trained, b_trained)
# test_accuracy = accuracy_score(y_test, y_pred_class)
# print(f"测试集准确率: {test_accuracy:.4f}")

可视化是理解模型如何工作的强大工具。我们可以绘制出数据的散点图，并叠加模型学习到的决策边界。

def plot_decision_boundary(X, y, w, b):
    """
    绘制二维特征下逻辑回归的决策边界
    """
    # 创建网格点
    x1_min, x1_max = X[:, 0].min() - 0.5, X[:, 0].max() + 0.5
    x2_min, x2_max = X[:, 1].min() - 0.5, X[:, 1].max() + 0.5
    xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max, 200),
                           np.linspace(x2_min, x2_max, 200))
    grid = np.c_[xx1.ravel(), xx2.ravel()]

    # 预测网格上每个点的概率
    _, Z_prob = predict(grid, w, b)
    Z = Z_prob.reshape(xx1.shape)

    # 绘制等高线和散点图
    plt.figure(figsize=(8, 6))
    plt.contourf(xx1, xx2, Z, levels=25, cmap='RdBu', alpha=0.6)
    plt.contour(xx1, xx2, Z, levels=[0.5], colors='black', linewidths=2)  # 决策边界
    plt.scatter(X[y.ravel()==0, 0], X[y.ravel()==0, 1], c='blue', label='Class 0', edgecolors='k')
    plt.scatter(X[y.ravel()==1, 0], X[y.ravel()==1, 1], c='red', label='Class 1', edgecolors='k')
    plt.xlabel('Feature 1')
    plt.ylabel('Feature 2')
    plt.title('Logistic Regression Decision Boundary')
    plt.legend()
    plt.colorbar(label='Predicted Probability of Class 1')
    plt.show()

# 使用一个简单的二维分类数据集进行演示
# from sklearn.datasets import make_classification
# X, y = make_classification(n_samples=100, n_features=2, n_informative=2, n_redundant=0, random_state=42)
# ... 训练模型 ...
# plot_decision_boundary(X, y, w_final, b_final)

4. 正则化：对抗过拟合的利器

当模型过于复杂（例如特征过多或多项式次数过高），它可能会完美“记忆”训练数据中的噪声和细节，导致在训练集上表现极好，但在未见过的测试集上表现糟糕。这种现象称为过拟合（或高方差）。正则化通过在损失函数中增加一个惩罚项，来约束模型参数的大小，从而降低模型复杂度，提高泛化能力。

4.1 L2正则化（岭回归）的原理与实现

最常用的正则化是L2正则化，它在原始损失函数的基础上，加上所有权重 w 的平方和（乘以一个正则化强度系数 lambda，通常写作 λ）。

线性回归的L2正则化损失函数： J(w, b) = (1/(2m)) * Σ(y_pred - y)^2 + (λ/(2m)) * Σ(w_j^2)

逻辑回归的L2正则化损失函数： J(w, b) = -(1/m) * Σ[y*log(f) + (1-y)*log(1-f)] + (λ/(2m)) * Σ(w_j^2)

注意，正则化项通常不包含偏置 b，因为 b 只是影响决策边界的位置，而不影响模型的弯曲复杂度。

实现正则化梯度下降时，只需要在权重 w 的梯度更新项中，额外加上 (λ/m) * w 这一项。

def compute_cost_linear_reg_with_regularization(X, y, w, b, lambda_):
    """带L2正则化的线性回归损失函数"""
    m = X.shape[0]
    y_pred = linear_hypothesis(X, w, b)
    mse_cost = compute_mse_cost(y, y_pred)
    # 正则化项，注意不惩罚b
    reg_cost = (lambda_ / (2 * m)) * np.sum(w ** 2)
    total_cost = mse_cost + reg_cost
    return total_cost

def gradient_descent_with_regularization(X, y, w, b, learning_rate, num_iterations, lambda_, model_type='linear'):
    """
    带L2正则化的梯度下降 (适用于线性和逻辑回归)
    model_type: 'linear' 或 'logistic'
    """
    m = X.shape[0]
    cost_history = []

    for i in range(num_iterations):
        if model_type == 'linear':
            y_pred = linear_hypothesis(X, w, b)
            error = y_pred - y
        elif model_type == 'logistic':
            z = np.dot(X, w) + b
            y_pred_prob = sigmoid(z)
            error = y_pred_prob - y
        else:
            raise ValueError("model_type must be 'linear' or 'logistic'")

        # 计算梯度
        dj_dw = (1 / m) * np.dot(X.T, error) + (lambda_ / m) * w  # 添加正则化项
        dj_db = (1 / m) * np.sum(error)  # b的梯度不变

        # 更新参数
        w = w - learning_rate * dj_dw
        b = b - learning_rate * dj_db

        # 计算并记录损失
        if model_type == 'linear':
            cost = compute_cost_linear_reg_with_regularization(X, y, w, b, lambda_)
        else:
            # 逻辑回归的带正则化损失计算
            z = np.dot(X, w) + b
            y_pred_prob = sigmoid(z)
            logistic_cost = compute_logistic_cost(y, y_pred_prob)
            reg_cost = (lambda_ / (2 * m)) * np.sum(w ** 2)
            cost = logistic_cost + reg_cost

        cost_history.append(cost)

    return w, b, cost_history

4.2 如何选择正则化参数 λ

λ 是一个超参数，需要手动调整。其选择是一个权衡：

• λ = 0：没有正则化，模型可能过拟合。
• λ 太大：惩罚过重，所有权重被压得非常接近0，模型会变得过于简单（欠拟合），可能变成一条水平线（对于线性回归）或一个常概率（对于逻辑回归）。

选择 λ 的黄金法则是使用验证集。将数据分为训练集、验证集和测试集。在训练集上用不同的 λ 值训练模型，然后在验证集上评估性能（例如，用准确率或F1分数），选择在验证集上表现最好的 λ。最后，用这个 λ 在完整的训练集上重新训练，并在独立的测试集上做最终评估。

def find_best_lambda(X_train, y_train, X_val, y_val, lambdas, model_type='logistic'):
    """
    通过验证集选择最佳的正则化参数 lambda
    """
    best_lambda = None
    best_val_accuracy = 0
    accuracies = []

    for lambda_ in lambdas:
        # 初始化参数
        w = np.random.randn(X_train.shape[1], 1) * 0.01
        b = 0.0

        # 使用当前lambda训练模型
        w_trained, b_trained, _ = gradient_descent_with_regularization(
            X_train, y_train, w, b, learning_rate=0.1, num_iterations=3000,
            lambda_=lambda_, model_type=model_type
        )

        # 在验证集上评估
        y_val_pred, _ = predict(X_val, w_trained, b_trained)
        val_accuracy = accuracy_score(y_val, y_val_pred)
        accuracies.append(val_accuracy)

        if val_accuracy > best_val_accuracy:
            best_val_accuracy = val_accuracy
            best_lambda = lambda_

    # 可视化不同lambda下的验证集准确率
    plt.plot(lambdas, accuracies, 'bo-')
    plt.xscale('log')  # lambda通常以对数尺度尝试
    plt.xlabel('Lambda (Regularization Strength)')
    plt.ylabel('Validation Accuracy')
    plt.title('Selecting Lambda via Validation Set')
    plt.grid(True)
    plt.show()

    print(f"最佳 lambda: {best_lambda}, 对应验证集准确率: {best_val_accuracy:.4f}")
    return best_lambda

# 示例用法
# lambdas_to_try = [0, 0.001, 0.01, 0.1, 1, 10, 100]
# best_lambda = find_best_lambda(X_train, y_train, X_val, y_val, lambdas_to_try, model_type='logistic')

在实际项目中，我经常发现 λ 在 [0.01, 0.1] 范围内效果不错，但这强烈依赖于数据和特征尺度。务必在标准化特征后再应用正则化，否则不同尺度的特征会受到不平等的惩罚。

从理论公式到可运行的代码，每一步都蕴含着对机器学习本质更深的理解。亲手实现这些基础算法，遇到的每一个bug，解决的每一个数值不稳定问题，都会让你对后续使用高级框架（如TensorFlow或PyTorch）时遇到的“黑箱”操作有更清晰的认知。当你能自如地调整学习率、观察损失曲线、通过正则化控制模型复杂度时，你才真正开始驾驭这些工具，而不仅仅是调用它们。