meshgrid
mesh[mɛʃ] 篩孔, 網線
numpy 的 meshgrid 是一個很好用的功能，可以產生每列 x 軸的座標及每行 y 軸的座標。
import numpy as np
a = range(5)
b = range(5)
x,y = np.meshgrid(a, b)
print("每列 x 的座標")
print(x)
print("每列 y 的座標")
print(y)
結果:
每列 x 的座標
[[0 1 2 3 4]
[0 1 2 3 4]
[0 1 2 3 4]
[0 1 2 3 4]
[0 1 2 3 4]]
每列 y 的座標
[[0 0 0 0 0]
[1 1 1 1 1]
[2 2 2 2 2]
[3 3 3 3 3]
[4 4 4 4 4]]
上述不論是 x 軸或 y 軸，都是[1,n] 的二維陣列。
標示座標
如果要標示每行每列的座標，代碼如下
import numpy as np
a = range(5)
b = range(5)
x,y = np.meshgrid(a, b)
for row in zip(x, y):
for col in zip(row[0], row[1]):
print(col, end='')
print()
結果:
(0, 0)(1, 0)(2, 0)(3, 0)(4, 0)
(0, 1)(1, 1)(2, 1)(3, 1)(4, 1)
(0, 2)(1, 2)(2, 2)(3, 2)(4, 2)
(0, 3)(1, 3)(2, 3)(3, 3)(4, 3)
(0, 4)(1, 4)(2, 4)(3, 4)(4, 4)
傳統思維產生格點
如果要在二維平面上產生格點如下圖
傳統思維的代碼如下
import pylab as plt
for x in range(20):
for y in range(20):
plt.scatter(x,y, s=0.5, c='b')
plt.show()
plt.scatter
plt.scatter(x, y) 其中的 x 是 x 軸的集合，y 是 y 軸的集合，代碼如下
import pylab as plt
x=[0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4]
y=[0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4]
plt.scatter(x,y,c='b')
plt.show()
上面藍色的 x, y 值，光 5*5 就打到快骨折，那麼如果是 20*20 的話，不就沒完沒了
meshgrid 產生 x , y 集合
把 meshgrid 產生的 x， y，由 [1,n] 的二維陣列轉成一維陣列 (x.reshape(-1))，再放入 plt.scatter(x, y) 就完成了。
import numpy as np import pylab as plt a = range(5) b = range(5) x, y = np.meshgrid(a, b) x = x.reshape(-1) y = y.reshape(-1) print(x) print(y) plt.scatter(x, y, s=0.5) plt.show() 結果: [0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4] [0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4]
plt 自動將二維轉成一維
上述的 x=x.reshape(-1) 其實是多餘的，因為 plt 會自動將多維扁平化成一維，所以只需如下代碼即可。
import numpy as np
import pylab as plt
a = range(5)
b = range(5)
x,y = np.meshgrid(a, b)
plt.scatter(x,y, s=0.5)
plt.show()
總結
meshgrid 的主要目的並不是要產生 (x, y) 的配對，而是要產生 x 軸的 list 及 y 軸的 list，方便傳給 plt.plot()
BGD 損失率圖解
BGD 是Bitch Gradient Descent (批次梯度下降) 的縮寫。
底下是產生迴歸線資料 Regression.py 的代碼
#pip install matplotlib
import numpy as np
np.random.seed(1)
def getData(n):
x = np.arange(-5, 5.1, 10 / n)
y = 3 * x + 2 + (np.random.rand(len(x)) - 0.5) * 20
return x, y
#取得等高線
def getContour(x, y):
a=np.arange(-10,16,1)
b=np.arange(-10,16,1)
mesh=np.meshgrid(a,b)
loss=np.zeros([len(a),len(b)])
for p_x, p_y in zip(x, y):
loss+=((mesh[0]*p_x+mesh[1])-p_y)**2
loss/=len(x)
return mesh,loss
BDG的損失函數梯度下降後，
$(\frac{\sigma Loss(a,b)}{\sigma a}=\frac{2}{n}\sum_{i=1}^{n}(ax_{i}+b-\tilde{y_{i}})*x_{i})$
$(\frac{\sigma Loss(a,b)}{\sigma b}=\frac{2}{n}\sum_{i=1}^{n}(ax_{i}+b-\tilde{y_{i}}))$
底下是BDG 類別 BDG.py 的代碼
import numpy as np
class BGD():
def __init__(self, a, b, x, y, lr):
self.a = a
self.b = b
self.x = x
self.y = y
self.lr = lr
self.a_old = a
self.b_old = b
self.loss = None;
def gradient(self):
#偏微分時，如果有10億筆資料，底下就要加總10億次，會很慢
#np.mean 跟 np.sum 一樣慢，把所有資料加總起來，再求平均數。
grad_a = 2 * np.mean((self.a * self.x + self.b - self.y) * (self.x))
grad_b = 2 * np.mean((self.a * self.x + self.b - self.y))
return grad_a, grad_b
def update(self):
# 計算梯度
grad_a, grad_b = self.gradient()
# 梯度更新
self.a_old = self.a
self.b_old = self.b
self.a = self.a - self.lr * grad_a
self.b = self.b - self.lr * grad_b
loss = ((self.a * self.x + self.b) - self.y) ** 2
self.loss = np.mean(loss)
底下是主程式的代碼
#BGD : Batch Gradient Descent : 批量梯度下降
import pylab as plt
from BGD import BGD
from Regression import *
np.random.seed(1)
epoch=50
x,y = getData(100)
mesh, contour=getContour(x,y)
fig,ax=plt.subplots(nrows=1, ncols=2, figsize=(12,4))
a2=ax[1].contourf(mesh[0], mesh[1], contour, 15, cmap=plt.cm.Purples)
ax[1].set_xlabel("a")
ax[1].set_ylabel("b")
plt.colorbar(a2, ax=ax[1])
lr = 0.1
init_a=-9#初始值
init_b=-9
ax[1].scatter(init_a, init_b, c='g')
gd=BGD(init_a, init_b, x, y, lr)
for i in range(epoch):
gd.update()
ax[0].clear()
ax[0].set_xlim(-5, 5)
ax[0].set_ylim(-30,30)
ax[0].scatter(x, y)
ax[0].plot([x[0], x[-1]], [gd.a*x[0]+gd.b, gd.a*x[-1]+gd.b], c="orange")
ax[0].set_title(f'{gd.a:.6f}x+{gd.b:.6f}')
ax[1].set_xlim(-10, 15)
ax[1].set_ylim(-10, 15)
ax[1].set_title(f'iter:{i+1:03d}, loss={gd.loss:.6f}')
ax[1].plot([gd.a_old, gd.a],[gd.b_old,gd.b],c='r')
ax[1].scatter(gd.a, gd.b, c='g')
plt.pause(0.01)
plt.show()
SGD 隨機梯度下降
SGD 為 Stochastic [stəˈkæstɪk] Gradient Descent 的縮寫。
BDG 損失函數梯度下降後，公式為$(\frac{\sigma Loss(a,b)}{\sigma a}=\frac{2}{n}\sum_{i=1}^{n}(ax_{i}+b-\tilde{y_{i}})*x_{i})$，必需將所有的 $(x_{i})$ 及 $(y_{i})$全都計算一遍。如果 i 只有幾百或幾千個，那電腦還可應付，但如果有上百萬，上億個點呢，這下就吃不消了。
SGD 就是把 $(\sum)$ 拿掉，然後使用亂數隨便取一個 x 及 y 來計算，這樣原本要加總計算上億個值，變成只要計算一次，就快很多了。當然啦，天下沒有白吃的午餐，如下圖所示，會呈現極不穩定的狀態。
SGD類別如下
from BGD import BGD
import numpy as np
class SGD(BGD):
def __init__(self, a, b, x, y, lr):
super().__init__(a, b, x, y, lr)
def gradient(self):
# 隨機取一筆資料進行計算
idx = np.random.randint(len(self.x))
grad_a = 2 * (self.a * self.x[idx] + self.b - self.y[idx]) * (self.x[idx])
grad_b = 2 * (self.a * self.x[idx] + self.b - self.y[idx])
return grad_a, grad_b
主程式如下
from Regression import *
import pylab as plt
from SGD import SGD
epoch=40
x,y=getData(100)
mesh, contour=getContour(x,y)
fig, ax=plt.subplots(nrows=1, ncols=2, figsize=(12,4))
a2=ax[1].contourf(mesh[0], mesh[1], contour,15, cmap=plt.cm.Purples)
plt.colorbar(a2,ax=ax[1])
lr = 0.01
init_a = -9; init_b = -9
ax[1].scatter(init_a, init_b, c='g')
gd = SGD(init_a, init_b, x, y, lr)
for i in range(epoch):
gd.update()
ax[0].clear()
ax[0].set_xlim(-5,5)
ax[0].set_ylim(-30, 30)
ax[0].scatter(x, y)
ax[0].plot([x[0], x[-1]],[gd.a*x[0]+gd.b, gd.a*x[-1]+gd.b], c="orange")
ax[0].set_title(f'{gd.a:.6f}x+{gd.b:.6f}')
print('iter=' + str(i) + ', loss=' + '{:.2f}'.format(gd.loss))
ax[1].set_xlim(-10,15)
ax[1].set_ylim(-10, 15)
ax[1].set_title(f'iter:{i:03d} Loss: {gd.loss:6f}')
ax[1].plot([gd.a_old, gd.a], [gd.b_old, gd.b], c='r')
ax[1].scatter(gd.a, gd.b, c='g')
ax[1].set_xlabel("a")
ax[1].set_ylabel("b")
plt.pause(0.01)
plt.show()
MBGD 小批量梯度下降
MBGD為 mini-batch gradient descent的縮寫。
因為SGD隨便選一個點來計算實在是不準確，所以就隨便選個100個點，或1000個點來計算。由底下的圖式可看出，確實平穩了許多，速度也快很多。MBGD 計算損失時，會產生 np.nan，所以如果遇到 np.nan 時，就中斷。
SGD 實在是上不了台面的，所以一般人家講的 SGD，其實是指 MBGD。
MBGD類別如下
from BGD import BGD import numpy as np class MBGD(BGD): def __init__(self, a, b, x, y, lr, batch_size): super().__init__(a, b, x, y, lr) self.batch_size=batch_size self.suffle=np.random.permutation(len(x)) self.refresh() self.update_batch() def refresh(self): self.suffle=np.random.permutation(len(self.x)) self.start=0 def update_batch(self): idx = self.suffle[self.start:self.start + self.batch_size] self.start+=self.batch_size self.x_batch = self.x[idx] self.y_batch = self.y[idx] def gradient(self): grad_a = 2 * np.mean((self.a * self.x_batch + self.b - self.y_batch) * (self.x_batch)) grad_b = 2 * np.mean((self.a * self.x_batch + self.b - self.y_batch))
#底下為提取批量到結尾，需重新打亂再重頭提取，避免 loss 為 nan
if self.start>len(self.x): self.refresh() self.update_batch() return grad_a, grad_b
主程式如下
from MBGD import MBGD
from Regression import *
import pylab as plt
epoch=200
x,y=getData(100)
mesh, contour=getContour(x,y)
fig, ax=plt.subplots(nrows=1, ncols=2, figsize=(12,4))
a2=ax[1].contourf(mesh[0], mesh[1], contour, 15, cmap=plt.cm.Purples)
plt.colorbar(a2,ax=ax[1])
lr=1e-2
init_a = -9; init_b = -9
ax[1].scatter(init_a, init_b, c='g')
batch_size=25
gd = MBGD(init_a, init_b, x, y, lr, batch_size)
for i in range(epoch):
gd.update()
ax[0].clear()
ax[0].set_xlim(-5,5)
ax[0].set_ylim(-30, 30)
ax[0].scatter(x, y)
ax[0].plot([x[0], x[-1]],[gd.a*x[0]+gd.b, gd.a*x[-1]+gd.b], c="orange")
ax[0].set_title(f'{gd.a:.6f}x+{gd.b:.6f}')
print('iter=' + str(i) + ', loss=' + '{:.2f}'.format(gd.loss))
ax[1].set_xlim(-10,15)
ax[1].set_ylim(-10, 15)
ax[1].set_title(f'iter:{i+1:03d} Loss: {gd.loss:6f}')
ax[1].plot([gd.a_old, gd.a], [gd.b_old, gd.b], c='r')
ax[1].scatter(gd.a, gd.b, c='g')
ax[1].set_xlabel("a")
ax[1].set_ylabel("b")
plt.pause(0.01)
plt.show()
SGD + 動量
上述的 MBGD 很容易卡在鞍點，所以可以再加入動量。本例應該是 MBGDM，但大家都稱 MBGD 為 SGD，所以本例命名為 SGDM
SGDM類別如下
from MBGD import MBGD import numpy as np class SGDM(MBGD): def __init__(self, a, b, x, y, lr, batch_size, gamma): super().__init__(a, b, x, y, lr, batch_size) self.gamma=gamma self.ma=0 self.mb=0 def update(self): grad_a, grad_b = self.gradient() self.a_old = self.a self.b_old = self.b self.ma = self.gamma * self.ma + self.lr * grad_a self.mb = self.gamma * self.mb + self.lr * grad_b self.a -=self.ma self.b -=self.mb loss = ((self.a * self.x + self.b) - self.y) ** 2 self.loss = np.mean(loss)
主程式如下
from SGDM import SGDM
from Regression import *
import pylab as plt
epoch=40
x,y=getData(100)
mesh, contour=getContour(x,y)
fig, ax=plt.subplots(nrows=1, ncols=2, figsize=(12,4))
a2=ax[1].contourf(mesh[0], mesh[1], contour, 15, cmap=plt.cm.Purples)
plt.colorbar(a2,ax=ax[1])
lr = 0.01
init_a = -9; init_b = -9
ax[1].scatter(init_a, init_b, c='g')
batch_size=5
gamma=0.8
gd = SGDM(init_a, init_b, x, y, lr, batch_size, gamma)
for i in range(epoch):
gd.update()
ax[0].clear()
ax[0].set_xlim(-5,5)
ax[0].set_ylim(-30, 30)
ax[0].scatter(x, y)
ax[0].plot([x[0], x[-1]], [gd.a * x[0] + gd.b, gd.a * x[-1] + gd.b], c="orange")
ax[0].set_title(f'{gd.a:.6f}x+{gd.b:.6f}')
print('iter=' + str(i) + ', loss=' + '{:.2f}'.format(gd.loss))
ax[1].set_xlim(-10,15)
ax[1].set_ylim(-10, 15)
ax[1].set_title(f'iter:{i+1:03d} Loss: {gd.loss:6f}')
ax[1].plot([gd.a_old, gd.a], [gd.b_old, gd.b], c= 'r')
ax[1].scatter(gd.a, gd.b, c='g')
ax[1].set_xlabel("a")
ax[1].set_ylabel("b")
plt.pause(0.01)
plt.show()
