人工智能學習筆記 python實作梯度下降法對多元函式求解

學習來自

梯度下降基本步驟如下圖所示

我們以一個二元函式為例計算

設一個二元函式為

y=0.5*(x1+x2)^2-x1*x2

一、則生成原函式影像代碼如下

#一、構建一個函式為 y=0.5*(x1+x2)^2-x1*x2的影像
#原函式如下
# 二維原始影像
def f2(x, y):
    return 0.15 * (x + 0.5) ** 2 + 0.25 * (y  - 0.25) ** 2 + 0.35 * (1.5 * x - 0.2 * y + 0.35 ) ** 2  

X1=np.arange(-4,4,0.2)
X2=np.arange(-4,4,0.2)
#Y = np.array(list(map(lambda t: f1(t),X)))
#Y = np.array(list(map(f2,zip(X1,X2))))
#print(Y)
X1, X2 = np.meshgrid(X1, X2) # 生成xv、yv，將X1、X2變成n*m的矩陣，方便后面繪圖
Y = np.array(list(map(f2,X1.flatten(),X2.flatten())))#這里壓縮成一維
Y.shape = X1.shape # 1600的Y圖還原成原來的（40,40）

生成的圖片如下

二、梯度下降步驟

1.隨機初始化引數值θ
2.計算梯度 這個點分別求關于x1、x2的偏導數，x1 =x1 - α*(dY/dx1),x2 =x2 - α*(dY/dx2)
3.修改引數值 alpha表示學習步長,也就是每次按照梯度減少的方向變化多少
4.按照3）迭代更新θ值，直至收斂或者θ 值的改變小于設定的閾值

# #初始化引數θ 閾值 和學習步長
x1 = 4
x2 = 4
alpha = 0.5
#保存梯度下降經過的點
GD_X1 = [x1]
GD_X2 = [x2]
GD_Y = [f2(x1,x2)]
# 定義y的變化量delta和迭代次數 閾值=limit
delta = 100 
iter_num = 0
limit = 1e-9 

# #計算梯度
def hx1(x, y):
    return 0.15 * 2 * (x + 0.5) + 0.25 * 2 * (1.5 * x - 0.2 * y + 0.35 ) * 1.5
def hx2(x, y):
    return 0.25 * 2 * (y  - 0.25) - 0.25 * 2 * (1.5 * x - 0.2 * y + 0.35 ) * 0.2

開始遞回

while delta > limit :
    tmp_x1 = x1 - alpha * hx1(x1,x2)
    tmp_x2 = x2 - alpha * hx2(x1,x2)
    tmp_y = f2(tmp_x1,tmp_x2)
    
    delta = np.abs(tmp_y - f2(x1,x2))
    x1 = tmp_x1
    x2 = tmp_x2
    GD_X1.append(x1)
    GD_X2.append(x2)
    GD_Y.append(tmp_y)
    iter_num += 1
    print(u"當前結果為:(%.5f, %.5f, %.5f)" % (x1, x2, tmp_y))
    print(u"迭代次數:%d" % iter_num)

迭代程序輸出記錄

三、畫圖步驟

fig = plt.figure(facecolor='w',figsize=(10,8))
ax = Axes3D(fig)
ax.plot_surface(X1,X2,Y,rstride=1,cstride=1,cmap=plt.cm.jet)
ax.plot(GD_X1,GD_X2,GD_Y,'ko-')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
 
ax.set_title(u'函式;\n學習率:%.3f; 最終解:(%.3f, %.3f, %.3f);迭代次數:%d' % (alpha, x1, x2, f2(x1,x2), iter_num))
plt.show()

得出最終影像為

梯度下降最終代碼

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
import random
from mpl_toolkits.mplot3d import Axes3D
import warnings
# 解決中文顯示問題
mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False

#一、構建一個函式為 y=0.5*(x1+x2)^2-x1*x2的影像
#原函式如下
# 二維原始影像
def f2(x, y):
    return 0.15 * (x + 0.5) ** 2 + 0.25 * (y  - 0.25) ** 2 + 0.35 * (1.5 * x - 0.2 * y + 0.35 ) ** 2  

X1=np.arange(-4,4,0.2)
X2=np.arange(-4,4,0.2)
#Y = np.array(list(map(lambda t: f1(t),X)))
#Y = np.array(list(map(f2,zip(X1,X2))))
#print(Y)
X1, X2 = np.meshgrid(X1, X2) # 生成xv、yv，將X1、X2變成n*m的矩陣，方便后面繪圖
Y = np.array(list(map(f2,X1,X2)))#這里壓縮成一維
Y.shape = X1.shape # 1600的Y圖還原成原來的（40,40）


# #二、梯度下降步驟
# # 1.隨機初始化引數值θ
# # 2.計算梯度 這個點分別求關于x1、x2的偏導數，x1 =x1 - α*(dY/dx1),x2 =x2 - α*(dY/dx2)  
# # 3.修改引數值  α表示學習步長,也就是每次按照梯度減少的方向變化多少
# # 4.按照3）迭代更新θ值，直至收斂或者θ 值的改變小于設定的閾值



# #計算梯度
def hx1(x, y):
    return 0.15 * 2 * (x + 0.5) + 0.25 * 2 * (1.5 * x - 0.2 * y + 0.35 ) * 1.5
def hx2(x, y):
    return 0.25 * 2 * (y  - 0.25) - 0.25 * 2 * (1.5 * x - 0.2 * y + 0.35 ) * 0.2

# #初始化引數θ 閾值 和學習步長
x1 = 4
x2 = 4
alpha = 0.5
#保存梯度下降經過的點
GD_X1 = [x1]
GD_X2 = [x2]
GD_Y = [f2(x1,x2)]
# 定義y的變化量和迭代次數
delta = f2(x1,x2)
iter_num = 0
limit = 1e-9 

while delta > limit :
    tmp_x1 = x1 - alpha * hx1(x1,x2)
    tmp_x2 = x2 - alpha * hx2(x1,x2)
    tmp_y = f2(tmp_x1,tmp_x2)
    
    delta = np.abs(tmp_y - f2(x1,x2))
    x1 = tmp_x1
    x2 = tmp_x2
    GD_X1.append(x1)
    GD_X2.append(x2)
    GD_Y.append(tmp_y)
    iter_num += 1
    print(u"當前結果為:(%.5f, %.5f, %.5f)" % (x1, x2, tmp_y))
    print(u"迭代次數:%d" % iter_num)

#print(GD_X1)
 
# 作圖
fig = plt.figure(facecolor='w',figsize=(10,8))
ax = Axes3D(fig)
ax.plot_surface(X1,X2,Y,rstride=1,cstride=1,cmap=plt.cm.jet)
ax.plot(GD_X1,GD_X2,GD_Y,'ko-')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
 
ax.set_title(u'函式;\n學習率:%.3f; 最終解:(%.3f, %.3f, %.3f);迭代次數:%d' % (alpha, x1, x2, f2(x1,x2), iter_num))
plt.show()