[吴恩达机器学习]1·梯度下降解一元线性回归

吴恩达机器学习系列课程：https://www.bilibili.com/video/BV164411b7dx

一元线性回归

给定数据集 \[ \left\{\left(x^{(i)},y^{(i)}\right),\;i=1,2,\cdots,m\right\} \] 试用一线性函数 \(y=\theta_1x+\theta_0\) 拟合之。当然，在数学上，我们熟知最小二乘法可以解决这个问题。

梯度下降法

这是我自学机器学习中遇到的第一个算法，其基本思想很简单。对于代价函数 \(J(\theta_0,\theta_1)\)，我们想得到它的一个极小值，只需要从任意点开始，选择函数 \(J\) 的梯度方向的逆方向，也即方向导数最大的方向移动一点点，然后不断重复这个过程。

我们知道梯度定义为： \[ \text{grad }J=\left\{\frac{\partial J}{\partial \theta_0},\frac{\partial J}{\partial \theta_1}\right\} \] 所以不断迭代进行： \[ \theta_j:=\theta_j-\alpha\cdot\frac{\partial J}{\partial\theta_j} \] 即可。其中，\(\alpha\) 就是这一步的长度，也称为学习率。

学习率的选取很重要，过小则梯度下降很慢，过大则有可能不收敛。

梯度下降解一元线性回归

应用在解一元线性回归问题上，我们定义代价函数为： \[ J(\theta_0,\theta_1)=\frac{1}{2m}\sum_{i=1}^m\left(\theta_1x^{(i)}+\theta_0-y^{(i)}\right)^2 \] 也即是每个数据纵坐标的真实值与估计值的差的平方和的平均，除以 \(2\) 仅是为了后续求导方便，没有什么本质的影响。

由于： \[ \text{grad }J=\left\{\frac{\partial J}{\partial \theta_0},\frac{\partial J}{\partial \theta_1}\right\}=\left\{\frac{1}{m}\sum_{i=1}^m\left(\theta_1x^{(i)}+\theta_0-y^{(i)}\right),\frac{1}{m}\sum_{i=1}^mx_i\left(\theta_1x^{(i)}+\theta_0-y^{(i)}\right)\right\} \] 于是每一次按照该方向的逆方向走一小步即可。

C++ 实现

对于这个问题，我们可以通过预处理达到每次迭代 \(O(1)\) 地更新。

#include<bits/stdc++.h>

using namespace std;

const double eps = 1e-8;
const double alpha = 0.01;

int m;
double sumx, sumy, sumxy, sumx2;
double th0, th1;

int main(){
	freopen("ex1data1.txt", "r", stdin);
	double x, y;
	while(scanf("%lf,%lf", &x, &y) != EOF)
		m++, sumx += x, sumy += y, sumx2 += x * x, sumxy += x * y;
	th0 = th1 = 0;
	while(1){
		double nth0, nth1;
		nth0 = th0 - alpha * 1.0 / m * (th1 * sumx - sumy + m * th0);
		nth1 = th1 - alpha * 1.0 / m * (th1 * sumx2 - sumxy + sumx * th0);
		if(fabs(th0 - nth0) < eps && fabs(th1 - nth1) < eps)	break;
		th0 = nth0, th1 = nth1;
	}
	printf("%f %f\n", th0, th1);
	return 0;
}

得到的结果是：\(\theta_0=-3.895775,\,\theta_1=1.193033\).

如果把 \(\alpha\) 设为 \(0.1\)，上面的程序无法得到结果，这是因为步长太大以至于无法收敛，除了人为调小 \(\alpha\) 以外，我尝试了让程序自适应地调整 \(\alpha\). 具体地，如果 \(\alpha\) 设置合理，我们的代价函数应该是一个随着迭代次数单调递减的函数，所以倘若代价函数出现了增加，我们就需要调小 \(\alpha\)，基于这个思想，我写下了如下的程序：

#include<bits/stdc++.h>

using namespace std;

const int N = 105;
const double eps = 1e-8;
double alpha = 100;

int m;
double x[N], y[N];
double th0, th1;

inline double gradJ(int k){
	double res = 0;
	for(int i = 1; i <= m; i++)
		res += (th1 * x[i] + th0 - y[i]) * (k == 1 ? x[i] : 1);
	return res / m;
}

inline double calc(){
	double res = 0;
	for(int i = 1; i <= m; i++)
		res += (th1 * x[i] + th0 - y[i]) * (th1 * x[i] + th0 - y[i]);
	return res / m / 2;
}

int main(){
	freopen("ex1data1.txt", "r", stdin);
	while(scanf("%lf,%lf", &x[0], &y[0]) != EOF)
		m++, x[m] = x[0], y[m] = y[0];
	th0 = th1 = 0;
	double preJ = 1e9, J = 0;
	while(1){
		double nth0, nth1;
		nth0 = th0 - alpha * gradJ(0);
		nth1 = th1 - alpha * gradJ(1);
		if(fabs(th0 - nth0) < eps && fabs(th1 - nth1) < eps)	break;
		th0 = nth0, th1 = nth1;
		J = calc();
		if(J > preJ)	alpha /= 5;
		preJ = J;
	}
	printf("%f\n%f %f\n", alpha, th0, th1);
	return 0;
}

最终得到的结果：\(\alpha=0.0064,\,\theta_0=-3.895772,\,\theta_1=1.193033\).

Python 实现

Python可以更方便地作图，也是现在机器学习中最热门的语言。

首先看一下原数据的散点图：

如下是学习率设置为 \(0.01\) 的 python 代码：

import numpy as np
import matplotlib.pyplot as plt

X = []
Y = []
m = 0
alpha = 0.01
eps = 1e-8
th0, th1 = 0, 0

def gradJ(k):
	res = 0
	for x, y in zip(X, Y):
		res += (th1 * x + th0 - y) * (x if k == 1 else 1)
	return res / m

def calc():
	res = 0
	for x, y in zip(X, Y):
		res += (th1 * x + th0 - y) * (th1 * x + th0 - y)
	return res / 2 / m

with open("ex1data1.txt", "r") as infile:
	data = infile.readlines()
	for line in data:
		x, y = line.split(',')
		X.append(float(x))
		Y.append(float(y))
		m += 1


while 1:
	nth0 = th0 - alpha * gradJ(0)
	nth1 = th1 - alpha * gradJ(1)
	if abs(th0 - nth0) < eps and abs(th1 - nth1) < eps:
		break
	th0 = nth0
	th1 = nth1
print(th0, th1)
plt.scatter(X, Y)
plt.plot([min(X), max(X)], [th0+min(X)*th1, th0+max(X)*th1], c="magenta")
plt.show()