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

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

Newton Raphson方法概述

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

实现步骤

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

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

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

其中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.01batch_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 @ firstreturn w, idef newton_curve(X, y, w):m = X.shape[0]h = h_theta(X, w)error = h - ycost = -(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)) @ Xreturn cost, second, firstdef 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方法的逻辑回归。每种方法都有其独特的应用场景和优势。

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

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

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

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

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