基于Python的机器学习系列（8）：Newton Raphson逻辑回归

        在本篇博文中，我们将探讨一种比传统梯度下降更高效的优化方法——Newton Raphson方法，并学习如何在逻辑回归中应用它。Newton Raphson方法通过利用二阶导数的曲率信息，快速地找到使代价函数最小化的参数。尽管这种方法在处理较小规模的数据集时速度极快，但随着特征数量的增加，它的性能可能会受到影响，因为需要计算矩阵的逆以及二阶导数。

Newton Raphson方法概述

        Newton Raphson方法是一种迭代算法，旨在通过对目标函数的二阶导数进行分析，找到函数的最小值点。与传统的梯度下降算法不同，Newton Raphson方法不仅考虑了一阶导数（斜率），还考虑了二阶导数（曲率），从而能够更快速地收敛到最优解。

实现步骤

        在本次实现中，我们将Newton Raphson方法应用于逻辑回归，具体步骤如下：

1. 数据准备

   • 添加截距项
   • 将输入矩阵 X、目标向量 y 和权重向量 w 进行正确的形状处理
     • X 的形状为 (m, n)，其中 m 是样本数量，n 是特征数量
     • y 的形状为 (m, )，为目标标签
     • w 的形状为 (n, )，为模型参数
   • 划分训练集和测试集
   • 特征缩放以加速收敛
   • 清理缺失数据（如果有）
   • （可选）特征工程

2. 预测和计算损失

                损失函数采用交叉熵损失：

                其中h 定义为sigmoid函数： 

        3.计算曲率方向

                使用如下公式计算曲率方向：

                其中 H 为二阶导数矩阵（即Hessian矩阵）。

        4.更新参数

                使用如下规则更新参数：

        5.迭代

                迭代执行步骤2-4，直到达到最大迭代次数 max_iter，或新旧损失之间的差异小于预定义的阈值。

数据准备

        首先，我们需要生成并处理数据，包括添加截距项、特征缩放等。以下代码展示了这一过程：

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=500, n_features=5, n_redundant=1, n_informative=4,
                             n_clusters_per_class=2, random_state=14)
plt.scatter(X[:, 0], X[:, 1], marker='o', c=y,
            s=25, edgecolor='k')

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)

# 特征缩放
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

# 添加截距项
intercept = np.ones((X_train.shape[0], 1))
X_train   = np.concatenate((intercept, X_train), axis=1)
intercept = np.ones((X_test.shape[0], 1))
X_test    = np.concatenate((intercept, X_test), axis=1)

Newton Raphson算法实现

        在实现Newton Raphson算法时，我们需要计算梯度和Hessian矩阵，并使用这些信息来更新权重向量。以下代码展示了这一过程：

def newton(X, y, max_iter=1000):
    w = np.zeros(X.shape[1])
    l_rate = 0.01
    batch_size = int(0.1 * X.shape[0])
    for i in range(max_iter):
        ix = np.random.randint(0, X.shape[0])
        batch_X = X[ix:ix+batch_size]
        batch_y = y[ix:ix+batch_size]
        cost, second, first = newton_curve(batch_X, batch_y, w)
        if i % 500 == 0:
            print(f"Cost at iteration {i}", cost)
        H_inverse = np.linalg.pinv(second)
        w = w - l_rate * H_inverse @ first
    return w, i

def newton_curve(X, y, w):
    m = X.shape[0]
    h = h_theta(X, w)
    error = h - y
    cost = -(np.sum(-y * np.log(h) - (1 - y) * np.log(1 - h)))
    first = (1/m) * np.dot(X.T, error)
    second = X.T @ np.diag((h) * (1-h)) @ X
    return cost, second, first

def sigmoid(x):
    return 1 / (1 + np.exp(-x))
        
def h_theta(X, w):
    return sigmoid(X @ w)

def output(pred):
    return np.round(pred)

w, i = newton(X_train, y_train, max_iter=5000)

计算准确率

        接下来，我们使用训练好的模型对测试集进行预测，并计算模型的准确率：

from sklearn.metrics import accuracy_score
yhat = output(h_theta(X_test, w))
accuracy_score(y_test, yhat)

结语

        在这一系列文章中，我们先介绍了二元逻辑回归，接着扩展到多元逻辑回归，最后讨论了Newton-Raphson方法的逻辑回归。每种方法都有其独特的应用场景和优势。

1. 二元逻辑回归：适用于只有两个类别的分类问题。通过对数几率（logit）函数的应用，我们可以将线性回归扩展到分类任务中。这个方法简单而有效，是大多数分类问题的基础。

2. 多元逻辑回归：当问题有多个类别时，多元逻辑回归则是首选。它通过使用softmax函数来处理多类别分类问题，是二元逻辑回归的自然扩展。

3. Newton-Raphson方法的逻辑回归：相比于传统的梯度下降法，Newton-Raphson方法提供了一种通过二阶导数（Hessian矩阵）来快速收敛的方法。虽然在特征数量较少时表现优异，但其在特征数量较多时由于需要计算Hessian矩阵的逆，可能会导致计算开销过大。

        每种方法在不同的场景下都有其独特的优势和局限性。了解并合理选择这些方法，可以帮助我们更好地解决分类问题。

        敬请期待下一篇博文：基于Python的机器学习系列（8）：朴素贝叶斯 - 高斯模型。

如果你觉得这篇博文对你有帮助，请点赞、收藏、关注我，并且可以打赏支持我！

欢迎关注我的后续博文，我将分享更多关于人工智能、自然语言处理和计算机视觉的精彩内容。

谢谢大家的支持！