正则化

课程与数据集

本文对应 Andrew Ng（吴恩达）《Machine Learning》正则化部分，是《逻辑回归》的延伸。代码使用的 ex2data2.txt（两次测试分数与芯片是否合格）为课程官方编程练习数据。

数据与参考代码可在黄海广教授整理的中文笔记仓库下载：fengdu78/Coursera-ML-AndrewNg-Notes（见 code/ex2/ 目录）。课程主页：Coursera Machine Learning。

---

一、为什么要进行正则化

1.1 过拟合

所谓过拟合，就是模型对数据拟合得很好，但是根本无法反映数据的变化规律。例如拉格朗日插值法（我强烈建议大家去 b 站看乐正垂心的那一期视频），它可以给出完美经过所有已知点的曲线，但这条曲线显然无法反映数据的变化趋势。一条可以经过所有数据点但是弯弯绕绕、次数和项过于复杂的曲线根本无法预测新的数据——也就是对原有数据拟合得好，但是没有任何预测价值。

1.2 欠拟合

欠拟合，顾名思义，就是曲线的复杂性根本满足不了拟合的需要。比如，一些明显呈现出对数变化规律的点，我们只会拿直线去拟合，这就达不到拟合源数据的需要了。

1.3 从线性到多项式

当我们将数据可视化之后，发现数据根本不能用最简单的线性回归或最简单的逻辑回归来拟合时，我们就要考虑将那个经典式子 $\overrightarrow{\theta }$ 从一次的线性式子，变成高次的多项式，例如将

$${\theta }_1{x}_1+{\theta }_2{x}_2+{\theta }_3{x}_3$$

转变为

$${\theta }_{11}{x}_1^2+{\theta }_{22}{x}_2^2+{\theta }_{33}{x}_3^2+{\theta }_{12}{x}_1{x}_2+{\theta }_{13}{x}_1{x}_3+{\theta }_{23}{x}_2{x}_3$$

这个拟合效果不就强了吗？但是我们又要考虑过拟合的问题，不是所有的项都有利，那么如何避免过拟合呢？方法就是正则化。

---

二、如何进行正则化

既然正则化是为了防止过拟合而存在的，那么就可以把所有项的系数都平方（防止正负抵消）都加入到代价函数里面：

$$J(\theta )=\frac{1}{m}\sum _{i=1}^m[-y^{(i)}\log (h_{\theta }(x^{(i)}))-(1-y^{(i)})\log (1-h_{\theta }(x^{(i)}))]+\frac{\lambda }{2m}\sum _{j=1}^n{\theta }_j^2$$

仔细观察这个式子，会发现如果系数过大，那么 cost 也会变大，这就是正则化的惩罚机制。它避免了各个项系数过大导致曲线过于弯曲的问题，让虽然不能适配每一个点，但是更有规律的那条曲线胜出。

---

三、代码

3.1 导入库与可视化

导入库和数据可视化不用多说：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.optimize as opt
from sklearn import linear_model

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

path = 'ex2data2.txt'
data2 = pd.read_csv(path, header=None, names=['Text1', 'Text2', 'Accepted'])

# 可视化数据
positive = data2.loc[data2['Accepted'].isin([1]), :]
negative = data2.loc[data2['Accepted'].isin([0]), :]

fig, ax = plt.subplots(figsize=(12, 8))
ax.scatter(positive['Text1'], positive['Text2'], s=50, c='b', marker='o', label='Accepted')
ax.scatter(negative['Text1'], negative['Text2'], s=50, c='r', marker='x', label='Rejected')
ax.legend()
ax.set_xlabel('Test 1 Score')
ax.set_ylabel('Test 2 Score')
plt.show()

3.2 构建多项式特征

我们来看看怎么构建多项式，此处选择最高次数为 5 次：

# 数据多项式拟合之前先做一步多项式处理
degree_max = 5
data2.insert(3, 'Ones', 1)
x1 = data2['Text1']
x2 = data2['Text2']

for i in range(1, degree_max):
    for j in range(0, i):
        data2['F' + str(i) + str(j)] = np.power(x1, i - j) * np.power(x2, j)  # 添加多项式的列
data2.drop('Text1', axis=1, inplace=True)  # 丢弃原来的一次列
data2.drop('Text2', axis=1, inplace=True)
print(data2.head())

# 运行完后，原来的数据表将会长这样

   Accepted  Ones       F10       F20       F21       F30       F31       F32       F40       F41       F42       F43
0         1     1  0.051267  0.002628  0.035864  0.000135  0.001839  0.025089  0.000007  0.000094  0.001286  0.017551
1         1     1 -0.092742  0.008601 -0.063523 -0.000798  0.005891 -0.043509  0.000074 -0.000546  0.004035 -0.029801
2         1     1 -0.213710  0.045672 -0.147941 -0.009761  0.031616 -0.102412  0.002086 -0.006757  0.021886 -0.070895
3         1     1 -0.375000  0.140625 -0.188321 -0.052734  0.070620 -0.094573  0.019775 -0.026483  0.035465 -0.047494
4         1     1 -0.513250  0.263426 -0.238990 -0.135203  0.122661 -0.111283  0.069393 -0.062956  0.057116 -0.051818

3.3 正则化代价函数

代价函数真心建议直接用 np.mean()，你拿 -np.sum(first + second) / len(x) 会算出 80 多：

def cost_reg(theta, x, y, learning_rate):
    first = y * np.log(sigmoid(x @ theta.T))
    second = (1 - y) * np.log(1 - sigmoid(x @ theta.T))
    reg = learning_rate / (2 * len(x)) * np.sum(np.power(theta[1:], 2))
    return -np.mean(first + second) + reg

# 数据准备
cols = data2.shape[1]
X = data2.iloc[:, 1:cols]
y = data2.iloc[:, 0:1]

X = np.array(X.values)
y = np.array(y.values)
theta_begin = np.zeros(11)

learning_rate = 1
cost1 = cost_reg(theta_begin, X, y, learning_rate)
print(cost1)

结果是：

0.6931471805599453

3.4 正则化梯度

你可以像上节课那样，计算第一个步长 grad1，只需要在函数里加上正则化项即可：

def gradient_reg(theta, x, y, learningRate):
    theta = theta.reshape(-1, 1)
    m = len(x)  # 样本数量
    error = sigmoid(x @ theta) - y  # 向量化计算误差
    # 计算梯度（向量化）
    grad = (x.T @ error) / m  # 基本梯度项
    # 添加 L2 正则化项（注意：不对 theta[0] 正则化）
    reg_term = (learningRate / m) * theta
    reg_term[0] = 0  # 跳过偏置项的正则化
    return (grad + reg_term).ravel()

grad1 = gradient_reg(theta_begin, X, y, learning_rate)
print(grad1)

结果：

[[0.00847458]
 [0.01878809]
 [0.05034464]
 [0.01150133]
 [0.01835599]
 [0.00732393]
 [0.00819244]
 [0.03934862]
 [0.00223924]
 [0.01286005]
 [0.00309594]]

3.5 优化求解

然后用：

result2 = opt.fmin_tnc(func=cost_reg, x0=theta_begin, fprime=gradient_reg, args=(X, y, learning_rate))
print(result2)

结果为：

(array([1.96323274e-03, 1.15183244e-03, -5.98648828e-03, -2.30810686e-03,
        4.69068291e-04, -9.02656339e-04, -1.64522767e-03, -4.53041378e-03,
        1.03585595e-05, -3.19517068e-03, -2.69567220e-04]), 92, 1)

3.6 预测准确率

再看看效果：

def predict(theta, x):
    probabilities = sigmoid(x @ theta.T)
    return (probabilities >= 0.5).astype(int)

theta_min = np.array(result2[0])
predictions = predict(theta_min, X)
accuracy = np.mean(predictions == y.flatten()) * 100
print(f'accuracy = {accuracy}%')

结果在 60% 左右就可以了。

3.7 可选：sklearn 实现

也可以用：

model = linear_model.LogisticRegression(penalty='l2', C=1.0)
r = model.fit(X, y.ravel())  # 一定要将 y 展平为一维
print(r)
print(model.score(X, y))

---

前置知识：《逻辑回归》笔记（ex2data1.txt）。数据集均来自 Coursera-ML-AndrewNg-Notes。