目 录
问题重述
附加问题
步骤实施
1.查看Scipy官网SciPy,找到优化有关的模块(Optimize)
2.研究多种优化策略,选择最符合代码的方案进行优化
3.minimize函数参数及其返回值
4.代码展示
5.结果展示
6.进一步优化
6.1对如下函数方法进行优化
6.2基准测试
6.3 发现
测试文件附录
任务清单
问题重述
在二维平面有n个点,如何画一条直线,使得所有点到该直线距离之和最短
如果能找到,请给出其损失函数
附加问题
1.使用Scipy优化上述问题
2.主代码中不得出现任何循环语法,出现一个扣10分
步骤实施
1.查看Scipy官网SciPy,找到优化有关的模块(Optimize)
2.研究多种优化策略,选择最符合代码的方案进行优化
名称 | 特点 | 应用场景 |
Scalar Functions Optimization | 用于最小化或最大化单个标量函数的,通常用于解决一维问题 | 目标函数只返回一个标量(单个值) |
Local (Multivariate) Optimization | 适用于多变量问题,需要梯度函数,不过会自动寻找梯度更新目标值 | 在参数空间中找到局部最小值或最大值 |
Global Optimization | 寻找函数的全局最小值或最大值,包含多个局部最值 | 在计算条件允许的条件下可以得到全局最优解 |
序号 | 名称 | 使用方法 | 适用条件 |
1 | Nelder-Mead | Nelder-Mead单纯形法 | 适用于一般的非线性问题 |
2 | Powell | Powell方法 | 适合多维非约束优化的方法 |
3 | CG | 共轭梯度法(Conjugate Gradient) | 适用于二次优化问题或大规模问题 |
4 | BFGS | 拟牛顿BFGS算法 | 适用于大多数非线性优化问题的常用方法,尤其是当梯度信息可用时 |
5 | Newton-CG | 牛顿共轭梯度法 | 适用于大多数非线性优化问题,但相对于BFGS需要更多的内存 |
6 | L-BFGS-B | 限制内存BFGS算法 | 适用于大规模问题,因为它限制了内存使用 |
7 | TNC | 截断牛顿法 | 适用于大多数非线性优化问题,并且能够处理约束条件 |
8 | COBYLA | 约束优化 | 适用于具有约束条件的问题 |
9 | SLSQP | 顺序最小二乘法 | 适用于具有约束条件的问题,并且能够处理线性和非线性约束 |
10 | trust-constr | 信任区域约束优化方法 | 适用于有约束条件的问题,并且可以处理线性和非线性约束 |
11 | dogleg | 信任域Dogleg方法 | 适用于具有约束条件的问题 |
12 | trust-ncg | 信任区域牛顿共轭梯度法 | 适用于约束优化问题 |
13 | trust-krylov | 信任区域Krylov子空间法 | 适用于约束优化问题 |
14 | trust-exact | 精确信任区域方法 | 适用于约束优化问题 |
此问题我们需要求最小值,所以我们采用minimize函数,并选择常用的BFGS策略
3.minimize函数参数及其返回值
原型如下:
scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)
挑五个主要的参数讲
1.fun:需要最小化的目标函数
这个函数应该接受一个输入向量,返回一个标量(单个值),表示损失函数的值。
2.x0:起始参数的初始猜测值
通常是一个数组或列表,表示参数的初始估计。
3.args:传递给目标函数的额外参数的元组
如果目标函数需要额外的参数,可以将它们作为元组传递给args参数。
4.method:选择优化方法的字符串
这是一个可选参数,如果未指定,默认使用'Nelder-Mead'方法。可以选择其他方法,
5.jac:表示目标函数的梯度(导数)的函数
如果提供了梯度函数,通常可以加速优化过程。如果不提供,优化算法会尝试数值估计梯度。
所以我们在优化代码的时候,
可以将 calcLoseFunction函数作为fun,
而k,b两个参数打成列表作为x0,
将XData,YData组成元组传递给arg
method选择BFGS
最后jac选择不写,便于对比两者速度差异
其返回值说明如下
1.x:优化的参数值。这是一个数组,包含找到的最优参数。
2.fun:最小化目标函数的最小值(损失函数的最小值)。
3.success:一个布尔值,表示优化是否成功收敛到最小值。
4.message:一个字符串,描述优化的终止消息。
5.nit:迭代次数,表示优化算法运行的迭代次数。
6.nfev:函数调用次数,表示评估目标函数的次数。
7.njev:梯度计算次数,表示计算目标函数梯度的次数(如果提供了梯度函数)。
8.hess_inv:Hessian矩阵的逆矩阵(如果提供了Hessian信息)。
9.jac:目标函数的梯度值。
4.代码展示
import numpy #发现直接用List就行了
import random
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from commonTools import *
# random.random()
# random.randint(start,stop)
#################全局数据定义区
# 数组大小
listSize=10
# 定义学习率 取尽量小0.001
learningRate=0.0001
#定义初始直线的 斜率k 和 截距b 45° 1单位距离
# 现在设置 k=0.5 检验程序
k,b=0.5,1
initialParams=[k,b]
#定义迭代次数
bfsNums=9999
#################全局数据定义区END
# 生成随机数
def generateRandomInteger(start, end):
# [1-100]
return random.randint(start, end)
# 打印本次随机生成的X,Y 便于快速粘贴复现
def printXYArray(XData,YData):
# 打印X
print("[", ",".join([str(i) for i in XData]), "]")
# 打印Y
print("[", ",".join([str(i) for i in YData]), "]")
#调用公共模块进行打印 便于快速查看粘贴
def printXYData(XData,YData):
loc=locals()
printArray(XData,loc)
printArray(YData,loc)
# 最小二乘法定义损失函数 并计算
#参考链接:https://blog.csdn.net/zy_505775013/article/details/88683460
# 求最小二乘法的最小值 最终结果应当是在learningRate一定情况下 这个最小的sum
def calcLoseFunction(params,XData,YData):
k, b = params
sum=0
for i in range(0,listSize):
# 使用偏离值的平方进行累和
sum+=(YData[i]-(k*XData[i]+b))**2
return sum
#梯度下降法
def calcGradientCorrection(b, k, XData, YData, learningRate, bfsNums):
for i in range(0, bfsNums):
sumk, sumb = 0, 0
for j in range(0, listSize):
# 定义预测值Y'
normalNum = k * XData[j] + b
# 计算逆梯度累和
sumk += -(2 / listSize) * (normalNum - YData[j]) * XData[j]
sumb += -(2 / listSize) * (normalNum - YData[j])
# 在逆梯度的方向上进行下一步搜索
k += learningRate * sumk
b += learningRate * sumb
return k, b
# 随机生成横坐标
XData=[generateRandomInteger(1,100) for i in range(listSize) ]
# 随机生成纵坐标
YData=[XData[i]+generateRandomInteger(-10,10) for i in range(listSize) ]
# 纯随机生成 但是可视化效果不直观
# YData=[generateRandomInteger(1,100) for i in range(listSize) ]
# 死值替换区
# XData=testArrayX
# YData=testArrayY
print("初始选取k={},b={}的情况下的损失函数值为sum={}".format(k,b,calcLoseFunction(initialParams,XData,YData)))
# 对k,b进行梯度修正
# k,b=calcGradientCorrection(b,k,XData,YData,learningRate,bfsNums)
#使用Scipy进行求解
result = minimize(calcLoseFunction, initialParams, args=(XData, YData), method='BFGS')
resultk,resultb=result.x
print("修正后:k={},b={},最小损失sum={},最小二乘法损失sums={}".format(resultk,resultb,result.fun,calcLoseFunction([resultk,resultb],XData,YData)))
print("调试数组")
printXYArray(XData,YData)
#画图
plt.plot(XData, YData, 'b.')
plt.plot(XData, resultk*numpy.array(XData)+resultb, 'r')
plt.show()
print("END")
5.结果展示
6.进一步优化
两个目标
1.优化损失函数中的for循环
2.对使用Scipy优化前后的代码进行基准测试,比较运行速度
6.1对如下函数方法进行优化
def calcLoseFunction(params,XData,YData):
k, b = params
sum=0
for i in range(0,listSize):
# 使用偏离值的平方进行累和
sum+=(YData[i]-(k*XData[i]+b))**2
return sum
使用numpy,优化后如下:
def calcLoseFunction(params,XData,YData):
XData,YData=np.array(XData),np.array(YData)
k, b = params
sum=np.sum((YData - (k * XData + b))**2)
return sum
无for优化后代码如下:
import numpy as np
import random
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from commonTools import *
#################全局数据定义区
# 数组大小
listSize=10
#定义初始直线的 斜率k 和 截距b 45° 1单位距离
# 现在设置 k=0.5 检验程序
k,b=0.5,1
initialParams=[k,b]
#################全局数据定义区END
# 生成随机数
def generateRandomInteger(start, end):
return random.randint(start, end)
#调用公共模块进行打印 便于快速查看粘贴
def printXYData(XData,YData):
loc=locals()
printArray(XData,loc)
printArray(YData,loc)
# 最小二乘法定义损失函数 并计算
def calcLoseFunction(params,XData,YData):
XData,YData=np.array(XData),np.array(YData)
k, b = params
sum=np.sum((YData - (k * XData + b))**2)
return sum
# 随机生成横坐标
XData=[generateRandomInteger(1,100) for i in range(listSize) ]
# 随机生成纵坐标
YData=[XData[i]+generateRandomInteger(-10,10) for i in range(listSize) ]
# 纯随机生成 但是可视化效果不直观
# YData=[generateRandomInteger(1,100) for i in range(listSize) ]
# 死值替换区
# XData=[ 49,74,62,54,20,14,27,74,23,50 ]
# YData=[ 47,65,56,57,21,21,32,81,27,46 ]
print("初始选取k={},b={}的情况下的损失函数值为sum={}".format(k,b,calcLoseFunction(initialParams,XData,YData)))
#使用Scipy进行求解
result = minimize(calcLoseFunction, initialParams, args=(XData, YData), method='BFGS')
resultk,resultb=result.x
print("修正后:k={},b={},最小损失sum={},最小二乘法损失sums={}".format(resultk,resultb,result.fun,calcLoseFunction([resultk,resultb],XData,YData)))
print("调试数组")
printXYData(XData,YData)
#画图
plt.plot(XData, YData, 'b.')
plt.plot(XData, resultk*np.array(XData)+resultb, 'r')
plt.show()
print("END")
其中公共模块commonTools.py 代码如下:
#########导包区
#########说明
#1.想要在公共模块区域使用变量列表 必须传进来 因为彼此的变量作用域不同
#########公共变量定义区
#这个locals应该是被引入的界面传进来,而不是从这拿
# loc=locals()
#########函数书写区
#1.获取变量名称
def getVariableName(variable,loc):
for k,v in loc.items():
if loc[k] is variable:
return k
#附带的打印变量名
def printValue(object,loc):
print("变量{}的值是{}".format(getVariableName(object,loc),object))
# 2.组装列表为字符串
def mergeInSign(dataList,sign):
# print(str(sign).join([str(i) for i in dataList]))
return str(sign).join([str(i) for i in dataList])
# 3.打印一个列表
def printArray(dataArray,loc):
print("列表{}的内容是:".format(getVariableName(dataArray,loc)),\
"[", ",".join([str(i) for i in dataArray]), "]"\
)
原先的代码如下:
import numpy #发现直接用List就行了
import random
import matplotlib.pyplot as plt
# random.random()
# random.randint(start,stop)
#################全局数据定义区
# 数组大小
listSize=10
# 定义学习率 取尽量小0.001
learningRate=0.0001
#定义初始直线的 斜率k 和 截距b 45° 1单位距离
# 现在设置 k=0.5 检验程序
k,b=0.5,1
#定义迭代次数
bfsNums=9999
#################全局数据定义区END
# 生成随机数
def generateRandomInteger(start, end):
# [1-100]
return random.randint(start, end)
# 打印本次随机生成的X,Y 便于快速粘贴复现
def printXYArray(XData,YData):
# 打印X
print("[", ",".join([str(i) for i in XData]), "]")
# 打印Y
print("[", ",".join([str(i) for i in YData]), "]")
# 最小二乘法定义损失函数 并计算
#参考链接:https://blog.csdn.net/zy_505775013/article/details/88683460
# 求最小二乘法的最小值 最终结果应当是在learningRate一定情况下 这个最小的sum
def calcLoseFunction(k,b,XData,YData):
sum=0
for i in range(0,listSize):
# 使用偏离值的平方进行累和
sum+=(YData[i]-(k*XData[i]+b))**2
return sum
#梯度下降法
def calcGradientCorrection(b, k, XData, YData, learningRate, bfsNums):
for i in range(0, bfsNums):
sumk, sumb = 0, 0
for j in range(0, listSize):
# 定义预测值Y'
normalNum = k * XData[j] + b
# 计算逆梯度累和 注意这里求偏导应当是两倍 不知道为什么写成1了
# 求MSE的偏导
sumk += -(2 / listSize) * (normalNum - YData[j]) * XData[j]
sumb += -(2 / listSize) * (normalNum - YData[j])
# 在逆梯度的方向上进行下一步搜索
k += learningRate * sumk
b += learningRate * sumb
return k, b
# 随机生成横坐标
XData=[generateRandomInteger(1,100) for i in range(listSize) ]
# 随机生成纵坐标
YData=[XData[i]+generateRandomInteger(-10,10) for i in range(listSize) ]
# 纯随机生成 但是可视化效果不直观
# YData=[generateRandomInteger(1,100) for i in range(listSize) ]
# 死值替换区
# XData=testArrayX
# YData=testArrayY
print("初始选取k={},b={}的情况下的损失函数值为sum={}".format(k,b,calcLoseFunction(k,b,XData,YData)))
# 对k,b进行梯度修正
k,b=calcGradientCorrection(b,k,XData,YData,learningRate,bfsNums)
print("修正后:k={},b={},最小损失sum={}".format(k,b,calcLoseFunction(k,b, XData, YData)))
print("调试数组")
printXYArray(XData,YData)
#画图
plt.plot(XData, YData, 'b.')
plt.plot(XData, k*numpy.array(XData)+b, 'r')
plt.show()
print("END")
到此,使用scipy并对for循环进行优化已经完成,下面我们使用程序对比优化后时间效率上有没有改进。
6.2基准测试
我们将先后代码的画图部分都注释
目录结构如下:
test.py代码如下:
import os #执行调用
import time #记录时间
DEBUG=False
execFileName="old.py" if DEBUG else "new.py"
if __name__=="__main__":
startTime = time.time()
os.system("python {}".format(execFileName))
endTime = time.time()
print("文件:{}执行耗时:{}ms".format(execFileName,endTime-startTime))
DEBUG为False
DEBUG为True
额,调用minimize函数在时间上不如自己写的梯度下降。。。。。
多次随机测试后发现结果依旧如此,可能是因为scipy引入了其他策略,导致了执行时间变长
6.3 发现
使用scipy在某种程度上可能能优化执行效率,但是在部分情况下可能耗时会略长于基本实现
测试文件附录
commonTools.py
#########导包区
#########说明
#1.想要在公共模块区域使用变量列表 必须传进来 因为彼此的变量作用域不同
#########公共变量定义区
#这个locals应该是被引入的界面传进来,而不是从这拿
# loc=locals()
#########函数书写区
#1.获取变量名称
def getVariableName(variable,loc):
for k,v in loc.items():
if loc[k] is variable:
return k
#附带的打印变量名
def printValue(object,loc):
print("变量{}的值是{}".format(getVariableName(object,loc),object))
# 2.组装列表为字符串
def mergeInSign(dataList,sign):
# print(str(sign).join([str(i) for i in dataList]))
return str(sign).join([str(i) for i in dataList])
# 3.打印一个列表
def printArray(dataArray,loc):
print("列表{}的内容是:".format(getVariableName(dataArray,loc)),\
"[", ",".join([str(i) for i in dataArray]), "]"\
)
test.py
import os #执行调用
import time #记录时间
DEBUG=True
execFileName="old.py" if DEBUG else "new.py"
if __name__=="__main__":
startTime = time.time()
os.system("python {}".format(execFileName))
endTime = time.time()
print("文件:{}执行耗时:{}ms".format(execFileName,endTime-startTime))
new.py
import numpy as np
import random
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from commonTools import *
#################全局数据定义区
# 数组大小
listSize=10
#定义初始直线的 斜率k 和 截距b 45° 1单位距离
# 现在设置 k=0.5 检验程序
k,b=0.5,1
initialParams=[k,b]
#################全局数据定义区END
# 生成随机数
def generateRandomInteger(start, end):
return random.randint(start, end)
#调用公共模块进行打印 便于快速查看粘贴
def printXYData(XData,YData):
loc=locals()
printArray(XData,loc)
printArray(YData,loc)
# 最小二乘法定义损失函数 并计算
def calcLoseFunction(params,XData,YData):
XData,YData=np.array(XData),np.array(YData)
k, b = params
sum=np.sum((YData - (k * XData + b))**2)
return sum
# # 随机生成横坐标
# XData=[generateRandomInteger(1,100) for i in range(listSize) ]
# # 随机生成纵坐标
# YData=[XData[i]+generateRandomInteger(-10,10) for i in range(listSize) ]
# 纯随机生成 但是可视化效果不直观
# YData=[generateRandomInteger(1,100) for i in range(listSize) ]
# 死值替换区
XData=[ 49,74,62,54,20,14,27,74,23,50 ]
YData=[ 47,65,56,57,21,21,32,81,27,46 ]
print("初始选取k={},b={}的情况下的损失函数值为sum={}".format(k,b,calcLoseFunction(initialParams,XData,YData)))
#使用Scipy进行求解
result = minimize(calcLoseFunction, initialParams, args=(XData, YData), method='BFGS')
resultk,resultb=result.x
print("修正后:k={},b={},最小损失sum={},最小二乘法损失sums={}".format(resultk,resultb,result.fun,calcLoseFunction([resultk,resultb],XData,YData)))
print("调试数组")
printXYData(XData,YData)
#画图
# plt.plot(XData, YData, 'b.')
# plt.plot(XData, resultk*np.array(XData)+resultb, 'r')
# plt.show()
print("END")
old.py
import numpy # 发现直接用List就行了
import random
import matplotlib.pyplot as plt
from commonTools import *
# random.random()
# random.randint(start,stop)
#################全局数据定义区
# 数组大小
listSize = 10
# 定义学习率 取尽量小0.001
learningRate = 0.0001
# 定义初始直线的 斜率k 和 截距b 45° 1单位距离
# 现在设置 k=0.5 检验程序
k, b = 0.5, 1
# 定义迭代次数
bfsNums = 9999
#################全局数据定义区END
# 生成随机数
def generateRandomInteger(start, end):
# [1-100]
return random.randint(start, end)
# 打印本次随机生成的X,Y 便于快速粘贴复现
def printXYArray(XData, YData):
# 打印X
print("[", ",".join([str(i) for i in XData]), "]")
# 打印Y
print("[", ",".join([str(i) for i in YData]), "]")
# 最小二乘法定义损失函数 并计算
# 参考链接:https://blog.csdn.net/zy_505775013/article/details/88683460
# 求最小二乘法的最小值 最终结果应当是在learningRate一定情况下 这个最小的sum
def calcLoseFunction(k, b, XData, YData):
sum = 0
for i in range(0, listSize):
# 使用偏离值的平方进行累和
sum += (YData[i] - (k * XData[i] + b)) ** 2
return sum
# 梯度下降法
def calcGradientCorrection(b, k, XData, YData, learningRate, bfsNums):
for i in range(0, bfsNums):
sumk, sumb = 0, 0
for j in range(0, listSize):
# 定义预测值Y'
normalNum = k * XData[j] + b
# 计算逆梯度累和 注意这里求偏导应当是两倍 不知道为什么写成1了
# 求MSE的偏导
sumk += -(2 / listSize) * (normalNum - YData[j]) * XData[j]
sumb += -(2 / listSize) * (normalNum - YData[j])
# 在逆梯度的方向上进行下一步搜索
k += learningRate * sumk
b += learningRate * sumb
return k, b
# 随机生成横坐标
XData = [generateRandomInteger(1, 100) for i in range(listSize)]
# 随机生成纵坐标
YData = [XData[i] + generateRandomInteger(-10, 10) for i in range(listSize)]
# 纯随机生成 但是可视化效果不直观
# YData=[generateRandomInteger(1,100) for i in range(listSize) ]
# 死值替换区
# XData=[ 49,74,62,54,20,14,27,74,23,50 ]
# YData=[ 47,65,56,57,21,21,32,81,27,46 ]
print("初始选取k={},b={}的情况下的损失函数值为sum={}".format(k, b, calcLoseFunction(k, b, XData, YData)))
# 对k,b进行梯度修正
k, b = calcGradientCorrection(b, k, XData, YData, learningRate, bfsNums)
print("修正后:k={},b={},最小损失sum={}".format(k, b, calcLoseFunction(k, b, XData, YData)))
print("调试数组")
printXYArray(XData, YData)
# 画图
# plt.plot(XData, YData, 'b.')
# plt.plot(XData, k * numpy.array(XData) + b, 'r')
# plt.show()
print("END")
任务清单
1.算法程序不使用任何for循环(已完成)
2.使用scipy对原先的代码进行优化(已完成)
3.对优化前后代码进行基准测试(已完成)