【机器学习练习 3】 - 一对多（One-vs-All, OvA）逻辑回归

本文链接：https://blog.youkuaiyun.com/wangshangshang09/article/details/146801711

【机器学习练习 3】 - 一对多（One-vs-All, OvA）逻辑回归

这段代码实现了一对多（One-vs-All, OvA）逻辑回归，用于多分类问题。

该代码涵盖了基于Python的解决方案，用于Coursera机器学习课程的第三个编程练习。有关详细说明和方程式。

数据集链接

对于此练习，我们将使用逻辑回归来识别手写数字（0到9）。我们将扩展我们在练习2中写的逻辑回归的实现，并将其应用于一对一的分类。让我们开始加载数据集。它是在MATLAB的本机格式，所以要加载它在Python，我们需要使用一个SciPy工具。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat

data = loadmat('ex3data1.mat')
data

data['X'].shape, data['y'].shape

((5000, 400), (5000, 1))

好的，我们已经加载了我们的数据。图像在martix X中表示为400维向量（其中有5,000个）。 400维“特征”是原始20 x 20图像中每个像素的灰度强度。类标签在向量y中作为表示图像中数字的数字类。

第一个任务是将我们的逻辑回归实现修改为完全向量化（即没有“for”循环）。这是因为向量化代码除了简洁外，还能够利用线性代数优化，并且通常比迭代代码快得多。但是，如果从练习2中看到我们的代价函数已经完全向量化实现了，所以我们可以在这里重复使用相同的实现。

sigmoid 函数

g 代表一个常用的逻辑函数（logistic function）为S形函数（Sigmoid function），公式为： $g\left( z \right)=\frac{1}{1+{{e}^{-z}}}$
合起来，我们得到逻辑回归模型的假设函数：
${{h}_{\theta }}\left( x \right)=\frac{1}{1+{{e}^{-{{\theta }^{T}}X}}}$

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

代价函数：
$J\left( \theta \right)=\frac{1}{m}\sum\limits_{i=1}^{m}{[-{{y}^{(i)}}\log \left( {{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)-\left( 1-{{y}^{(i)}} \right)\log \left( 1-{{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)]}$

def cost(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    reg = (learningRate / (2 * len(X))) * np.sum(np.power(theta[:,1:theta.shape[1]], 2))
    return np.sum(first - second) / len(X) + reg

如果我们要使用梯度下降法令这个代价函数最小化，因为我们未对 ${{\theta }_{0}}$ 进行正则化，所以梯度下降算法将分两种情形：
$\begin{align} & Repeat\text{ }until\text{ }convergence\text{ }\!\!\{\!\!\text{ } \\ & \text{ }{{\theta }_{0}}:={{\theta }_{0}}-a\frac{1}{m}\sum\limits_{i=1}^{m}{[{{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}}]x_{_{0}}^{(i)}} \\ & \text{ }{{\theta }_{j}}:={{\theta }_{j}}-a\frac{1}{m}\sum\limits_{i=1}^{m}{[{{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}}]x_{j}^{(i)}}+\frac{\lambda }{m}{{\theta }_{j}} \\ & \text{ }\!\!\}\!\!\text{ } \\ & Repeat \\ \end{align}$

以下是原始代码是使用for循环的梯度函数：

def gradient_with_loop(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    
    error = sigmoid(X * theta.T) - y
    
    for i in range(parameters):
        term = np.multiply(error, X[:,i])
        
        if (i == 0):
            grad[i] = np.sum(term) / len(X)
        else:
            grad[i] = (np.sum(term) / len(X)) + ((learningRate / len(X)) * theta[:,i])
    
    return grad

向量化的梯度函数

def gradient(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    parameters = int(theta.ravel().shape[1])
    error = sigmoid(X * theta.T) - y
    
    grad = ((X.T * error) / len(X)).T + ((learningRate / len(X)) * theta)
    
    # intercept gradient is not regularized
    grad[0, 0] = np.sum(np.multiply(error, X[:,0])) / len(X)
    
    return np.array(grad).ravel()

现在我们已经定义了代价函数和梯度函数，现在是构建分类器的时候了。对于这个任务，我们有10个可能的类，并且由于逻辑回归只能一次在2个类之间进行分类，我们需要多类分类的策略。在本练习中，我们的任务是实现一对一全分类方法，其中具有k个不同类的标签就有k个分类器，每个分类器在“类别 i”和“不是 i”之间决定。我们将把分类器训练包含在一个函数中，该函数计算10个分类器中的每个分类器的最终权重，并将权重返回为k X（n + 1）数组，其中n是参数数量。

from scipy.optimize import minimize

def one_vs_all(X, y, num_labels, learning_rate):
    rows = X.shape[0]
    params = X.shape[1]
    
    # k X (n + 1) array for the parameters of each of the k classifiers
    all_theta = np.zeros((num_labels, params + 1))
    
    # 添加偏置项
    X = np.insert(X, 0, values=np.ones(rows), axis=1) 
    
    # labels are 1-indexed instead of 0-indexed
    for i in range(1, num_labels + 1):
        theta = np.zeros(params + 1)
        y_i = np.array([1 if label == i else 0 for label in y])
        y_i = np.reshape(y_i, (rows, 1))
        
        # minimize the objective function
        # minimize 是 scipy.optimize 提供的最优化函数，使用 共轭梯度法（TNC） 进行优化。
        fmin = minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient)
        all_theta[i-1,:] = fmin.x
    
    return all_theta

这里需要注意的几点：首先，我们为theta添加了一个额外的参数（与训练数据一列），以计算截距项（常数项）。其次，我们将y从类标签转换为每个分类器的二进制值（要么是类i，要么不是类i）。最后，我们使用SciPy的较新优化API来最小化每个分类器的代价函数。如果指定的话，API将采用目标函数，初始参数集，优化方法和jacobian（渐变）函数。然后将优化程序找到的参数分配给参数数组。

实现向量化代码的一个更具挑战性的部分是正确地写入所有的矩阵，保证维度正确。

rows = data['X'].shape[0]
params = data['X'].shape[1]

all_theta = np.zeros((10, params + 1))

X = np.insert(data['X'], 0, values=np.ones(rows), axis=1)

theta = np.zeros(params + 1)

y_0 = np.array([1 if label == 0 else 0 for label in data['y']])
y_0 = np.reshape(y_0, (rows, 1))

X.shape, y_0.shape, theta.shape, all_theta.shape

((5000, 401), (5000, 1), (401,), (10, 401))

注意，theta是一维数组，因此当它被转换为计算梯度的代码中的矩阵时，它变为（1×401）矩阵。我们还检查y中的类标签，以确保它们看起来像我们想象的一致。

np.unique(data['y'])#看下有几类标签

array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10], dtype=uint8)

让我们确保我们的训练函数正确运行，并且得到合理的输出。

all_theta = one_vs_all(data['X'], data['y'], 10, 1)
all_theta

我们现在准备好最后一步 - 使用训练完毕的分类器预测每个图像的标签。对于这一步，我们将计算每个类的类概率，对于每个训练样本（使用当然的向量化代码），并将输出类标签为具有最高概率的类。

def predict_all(X, all_theta):
    rows = X.shape[0]
    params = X.shape[1]
    num_labels = all_theta.shape[0]
    
    # same as before, insert ones to match the shape
    X = np.insert(X, 0, values=np.ones(rows), axis=1)
    
    # convert to matrices
    X = np.matrix(X)
    all_theta = np.matrix(all_theta)
    
    # compute the class probability for each class on each training instance
    h = sigmoid(X * all_theta.T)
    
    # create array of the index with the maximum probability
    h_argmax = np.argmax(h, axis=1)
    
    # because our array was zero-indexed we need to add one for the true label prediction
    h_argmax = h_argmax + 1
    
    return h_argmax

现在我们可以使用predict_all函数为每个实例生成类预测，看看我们的分类器是如何工作的。

y_pred = predict_all(data['X'], all_theta)
correct = [1 if a == b else 0 for (a, b) in zip(y_pred, data['y'])]
accuracy = (sum(map(int, correct)) / float(len(correct)))
print ('accuracy = {0}%'.format(accuracy * 100))

accuracy = 94.46%

note

训练多个分类器的步骤是 One-vs-All 方法的核心，它的关键思想是：

将 多分类问题（num_labels 类）转化为 多个二分类问题。
对每个类别 ( i ) 训练一个分类器：
- 让当前类别 ( i ) 的标签变成 正类（1），所有其他类别变成 负类（0）。
- 训练一个 逻辑回归模型 来区分当前类别与其他类别。
- 依次重复这个过程，共训练 num_labels 个分类器。

代码逐行解析

for i in range(1, num_labels + 1):

循环遍历每个类别（类别标签 1, 2, ..., num_labels）。
这里 i 代表当前正在训练的类别编号。

theta = np.zeros(params + 1)

初始化参数 theta，形状为 (n+1, )，对应于逻辑回归的权重，包括 theta_0（偏置项）。

y_i = np.array([1 if label == i else 0 for label in y])
y_i = np.reshape(y_i, (rows, 1))

y_i 是一个二分类标签：
- 如果 y == i，设为 1（当前类别）。
- 否则设为 0（不是当前类别）。

这样，我们把 y 转换成了一个二分类问题，例如：

原始 y = [1, 2, 1, 3, 2]
训练 i=1 时: y_i = [1, 0, 1, 0, 0]
训练 i=2 时: y_i = [0, 1, 0, 0, 1]
训练 i=3 时: y_i = [0, 0, 0, 1, 0]

这样，我们针对每个类别 ( i ) 分别训练一个二分类器。

fmin = minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient)

minimize 函数用于优化逻辑回归的损失函数，找到最优 theta。
传入 cost（损失函数），gradient（梯度），以及 X, y_i, learning_rate 作为参数。

all_theta[i-1,:] = fmin.x

fmin.x 是优化后得到的 theta（最佳权重）。
存入 all_theta[i-1, :]，对应当前类别 ( i ) 的分类器参数。

示例

假设我们有 3 类数据，输入数据 X 为：

X = [[0.1, 0.2], 
     [0.4, 0.5], 
     [0.7, 0.8], 
     [1.0, 1.1]]

标签 y：

y = [1, 2, 3, 1]

模型会进行 3 次循环：

训练 i=1 的分类器（类别 1 vs 其他）：
- y_i = [1, 0, 0, 1]
- 训练逻辑回归，得到 theta_1。
训练 i=2 的分类器（类别 2 vs 其他）：
- y_i = [0, 1, 0, 0]
- 训练逻辑回归，得到 theta_2。
训练 i=3 的分类器（类别 3 vs 其他）：
- y_i = [0, 0, 1, 0]
- 训练逻辑回归，得到 theta_3。

最终得到：

all_theta = [theta_1,   # 1类的参数
             theta_2,   # 2类的参数
             theta_3]   # 3类的参数

用于分类时，每个 X 都会计算 theta_i 的预测概率，并选取最大值的类别作为最终预测。

总结

为什么要训练多个分类器？

逻辑回归本身是二分类模型，不能直接用于多分类。
One-vs-All 方法：每个类别训练一个二分类模型。
训练后，我们有 num_labels 个分类器，每个针对一个类别。

代码核心逻辑

遍历每个类别 i：
- 让 y_i 变成二分类（1=当前类别, 0=其他类别）。
- 训练二分类逻辑回归模型，求解 theta。
- 将 theta 存入 all_theta。
最终得到 all_theta，可用于预测新数据的类别。

这样，每个类别都有一个专属分类器，最终可用于分类预测！🚀