实验题目
- Generate n = 2,000 points uniformly at random in the two-dimensional unit square. Which point do you expect the centroid to be?
- What objective does the centroid of the points optimize?
- Apply gradient descent (GD) to find the centroid.
- Apply stochastic gradient descent (SGD) to find the centroid. Can you
say in simple words, what the algorithm is doing?
实验要求
1)编程语言不限。
2)作业包含一份报告(word 或 pdf 格式)及代码加注释，并打包到.zip，其中 zip 文件的命名格式为学号_姓名。
3)不允许使用梯度下降相关的库函数。
4)禁止抄袭。
实验过程及代码
Q1
Generate n = 2,000 points uniformly at random in the two-dimensional unit square. Which point do you expect the centroid to be?
在数学和物理学中，平面图形的质心或几何中心是图中所有点的算术平均值；在几何学中，重心时质心的同义词；在天体物理学和天文学中，重心是两个或多个轨道相互环绕的物质的质心。如果物理对象具有均匀的密度，则其质心与其形状的质心相同。
首先生成2000个分布在0~1范围之间的随机散点，从统计学的角度理解，它们具有均匀的密度，质心是最能代表这群点的点，它到所有点距离之和最小。由于点的生成是随机的，所以期望的点是(0.5, 0.5)。
我们利用如下的代码生成2000个分布在0~1范围之间的随机散点，用蓝色点标出，同时用红点标注我们期望的质心。
import numpy as np
import matplotlib.pyplot as plt
Sample_num = 2000
x = np.random.random_sample((Sample_num, ))
y = np.random.random_sample((Sample_num, ))
plt.figure(figsize=(8, 8))
plt.plot(x, y, 'b.')
plt.plot(0.5, 0.5, 'ro')
Q2
What objective does the centroid of the points optimize?
在*Q1*中我们分析得到，质心是一群点中到其他所有点距离和最小的点，所以我们的目标就是最小化距离之和。
假设当前点为:C(x,y)，与其他所有点的距离为dist(x,y)。所以优化的目标是：
此时我们可以定义损失函数如下：
def cost(c, points):
'''
params:
c: (x,y) 表示待求的质心坐标
points: 所有的点
return:
dist: c到所有点的距离之和
'''
return sum(sum((c - points) ** 2, axis=1) ** 0.5)Q3
Apply gradient descent (GD) to find the centroid.
损失函数
根据损失函数cost，求出其在当前点(x, y)两个维度上的梯度(dx, dy)，梯度求解的公式如下：
其代码如下：
def gradient_raw(c, points):
'''
params:
c: (x,y) 表示待求的质心坐标
points: 所有的点
return:
(dx, dy), 在c点时cost在两个维度方向的梯度
'''
x = points[:, 0]
y = points[:, 1]
dist = sum((c - points) ** 2, axis=1) ** 0.5
dx = sum((c[0] - x) / dist) #求x偏导数
dy = sum((c[1] - y) / dist) #求y偏导数
return np.asarray([dx, dy])GD
这一步设置一些参数，编写好梯度下降函数。
def gradientDescent_raw(points, start=array([1, 1]), theta=0.01, iterations=500, eps=1e-6):
'''
params:
start: 初始化点，即出发点，为了便于观察，可以设置为四个顶点中的任意一个，这里默认为右上角的顶点
points: 所有采样点
theta: 步长，即学习率，控制步长
iterations: 迭代次数
eps: 阈值，通过比较前后连续两次的cost差值控制迭代。
return:
x_route: centroid路径
x: 收敛时的centroid
i: end iteration, 退出迭代时已经迭代了多少次
'''
x = start
x_route = x
pre_cost = cost(x, points) # 初始cost
for i in range(iterations):
grad = gradient(x, points)
x_i = x - theta * grad
cur_cost = cost(x_i, points)
if abs(pre_cost - cur_cost) > eps: # 前一个cost与当前cost的差大于阈值，可以继续迭代
x = x_i
pre_cost = cur_cost
else:
return x_route, x, i+1
x_route = vstack((x_route, x))
return x_route, x, i+1结果绘图与分析
编写辅助函数与测试函数。
def plotRoute(points, route, c, endIteration):
plt.plot(points[:, 0], points[:, 1], 'b,')
plt.plot(route[:, 0], route[:, 1], 'r.')
plt.plot(route[:, 0], route[:, 1], 'k-')
plt.xlabel('c = (%.3f, %.3f), endAt = %d' % (c[0], c[1], endIteration))
print(c)
plt.show()
def test():
x_b, c, i = gradientDescent(points, theta=0.1,, iterations=1000, lr_decay=0.8)
plotRoute(points, x_b, c, i)运行上面的程序并绘图。
结果分析：
问题1：不收敛
图示路径与最终得到的centroid(x,y)表明达到最大迭代次数后，没有收敛，尝试了增大迭代次数到10000，100000之后依然得到的是类似的结果。
观察发现路径不是一点一点慢慢向中心靠拢，而是在两边很夸张地跳跃，因此问题出在theta，也就是我们熟悉的学习率上。当某一步“跨”的过大，即本次迭代更新后的cur_cost比pre_cost要大的时候，说明我们不是朝着下降的方向在走，此时要调整我们的步伐，迈小一点的步子。因此我们更新我们的梯度下降函数如下。
def gradientDescent(points, start=array([1, 1]), theta=0.01, iterations=500, eps=1e-6, lr_decay=0.5):
'''
params:
start: 初始化点，即出发点，为了便于观察，可以设置为四个顶点中的任意一个，这里默认为右上角的顶点
points: 所有采样点
theta: 步长，即学习率，控制步长
iterations: 迭代次数
eps: 阈值，通过比较前后连续两次的cost差值控制迭代
lr_decay: 步长衰减比例
return:
x_route: centroid路径
x: 收敛时的centroid
i: end iteration, 退出迭代时已经迭代了多少次
'''
x = start
x_route = x
pre_cost = cost(x, points) # 初始cost
for i in range(iterations):
grad = gradient(x, points)
x_i = x - theta * grad
cur_cost = cost(x_i, points)
if pre_cost - cur_cost > eps: # 前一个cost与当前cost的差大于阈值，可以继续迭代
x = x_i
pre_cost = cur_cost
elif cur_cost - pre_cost > eps: # 这一步“跨”的过大，调小步长
theta = theta * lr_decay
else:
return x_route, x, i+1
x_route = vstack((x_route, x))
return x_route, x, i+1再次运行结果绘图。
出现新的问题。
问题2:第一步直接“跨”到中心位置附近
在gradient_raw函数中，我们得到的梯度不仅包含方向，也包含了在该方向上的长度， 这个长度数值可能较大，因此我们将其单位化，即dx和dy分别除以它们的欧式距离。因此更新我们的gradient函数如下。
def gradient(c, points):
'''
params:
c: (x,y) 表示待求的质心坐标
points: 所有的点
return:
(dx, dy), 在c点时cost在两个维度方向的梯度
'''
x = points[:, 0]
y = points[:, 1]
dist = sum((c - points) ** 2, axis=1) ** 0.5
dx = sum((c[0] - x) / dist) #求x偏导数
dy = sum((c[1] - y) / dist) #求y偏导数
s = (dx ** 2 + dy ** 2) ** 0.5 # 欧式距离
dx = dx/s
dy = dy/s
return np.asarray([dx, dy])运行结果绘图。
图中红色的点(由于比较密集，看起来像线)，表示路径上的点，在稳步向中心靠拢，最终得到的质心也接近理想的质心。至此，完成本部分实验。
Q4
Apply stochastic gradient descent (SGD) to find the centroid. Can you say in simple words, what the algorithm is doing?
随机梯度下降和梯度下降的区别：
1. 梯度下降需要所有点都参与梯度更新；随机梯度下降是每次迭代时，随机选取一个点计算梯度。
2. 梯度下降最终的结果是一般是稳定的，即稳定收敛到某个点；随机梯度下降可能会陷入某个局部最小值。
修改gradient参数，使得迭代时随机选取一个点。
def stochasticgradient(c, points):
'''
params:
c: (x,y) 表示待求的质心坐标
points: 所有的点
return:
(dx, dy), 在c点时cost在两个维度方向的梯度
'''
r = choice(points)
dx = (c[0] - r[0]) / sum((c - r) ** 2) ** 0.5 # x偏导
dy = (c[1] - r[1]) / sum((c - r) ** 2) ** 0.5 # y偏导
s = (dx ** 2 + dy ** 2) ** 0.5
dx = dx/s
dy = dy/s
return array([dx, dy])我们修改gradientDescent函数，在计算梯度的时候调用stochasticgradient方法，这样我们把计算梯度所用的函数作为一个参数，同样传递进去。修改后代码如下：
def gd_final(points, gd_method=gradient,start=array([1, 1]), theta=0.01, iterations=500, eps=1e-6, lr_decay=0.5):
'''
params:
gd_method: 迭代的方式，可以在(stochasticgradient, gradient)中选择
start: 初始化点，即出发点，为了便于观察，可以设置为四个顶点中的任意一个，这里默认为右上角的顶点
points: 所有采样点
theta: 步长，即学习率，控制步长
iterations: 迭代次数
eps: 阈值，通过比较前后连续两次的cost差值控制循环结束。
return:
x_route: centroid路径
x: 收敛时的centroid
i: end iteration, 退出迭代时已经迭代了多少次
'''
x = start
x_route = x
pre_cost = cost(x, points) # 初始cost
for i in range(iterations):
grad = gd_method(x, points)
x_i = x - theta * grad
cur_cost = cost(x_i, points)
if pre_cost - cur_cost > eps: # 前一个cost与当前cost的差大于阈值，可以继续迭代
x = x_i
pre_cost = cur_cost
elif cur_cost - pre_cost > eps: # 这一步“跨”的过大，调小步长
theta = theta * lr_decay
else:
return x_route, x, i+1
x_route = vstack((x_route, x))
return x_route, x, i+1最终得到下图：
这张图是在多次调整参数多次运行之后得到结果相对较好的图。 下面两张则是陷入局部最小时的情况。
可以看到，由于是随机选取的点，所以在找最优值的过程是曲折的。但是最终还是会收敛。