如何在python中加速非线性优化的多次迭代？
How to speed up the multiple iterations of non-linear optimization in python?
我正在尝试使用 PyGMO package 求解非线性优化
优化 class 是单独定义的，然后通过单独的函数 dyn_optimGMO 调用。必须为变量（节点）定义的 1000 个随机初始向量完成并保存此优化 inits ( or init_val)
使用 timeit 模块我发现每次迭代大约需要 17 seconds 才能完成。这意味着 1000 iterations 大约需要 5 hours 。这是一个非常大的时间。
如果我必须对 20 perturb 个节点重复此操作，那么总迭代次数将达到 200000，这将花费上面计算的线性时间。
我尝试通过使用 python multiprocessing 模块为 20 个扰动节点中的每一个并行化每组 1000 次迭代来解决这个问题。但这于事无补。
我也尝试过使用 Numba jit 函数，但它们无法识别 pyGMO 模块，因此失败了。
有什么方法可以并行化此代码并使其在任意次数的迭代中都更快？
如果我的问题足够清楚，请告诉我，如果不够清楚，我会根据需要添加详细信息。
import numpy as np
import pygmo as pg
matL = np.random.rand(300,300) ; node_len = 300
inits = []; results = []
perturb = {25:0} #setting a random node, say, node 25 to 0
class my_constrained_udp:
def __init__(self):
pass
def fitness(self, x):
matA = np.matrix(x)
obj1 = matA.dot(matL).dot(matA.T)[0,0] #this is int value
ce1 = sum(init_val) - sum(x)
return [obj1, ce1]
def n_objs(self): # no of objectives
return 1
def get_nec(self): #no of equality constraints
return 1
def get_nic(self): #no of in-equality constraints
return 0
def get_bounds(self): #lower and upper bounds: use this to perturb the node
lowerB = np.array([0]*node_len); upperB = np.array([1]*node_len)
if perturb:
for k,v in perturb.items():
lowerB[k] = v; upperB[k] = v
return (lowerB,upperB)
def gradient(self, x):
return pg.estimate_gradient_h(lambda x: self.fitness(x), x)
def dyn_optimGMO(matL, node_len ,init):
if perturb:
for k,v in perturb.items(): init_val[k] = v #setting perturbations in initial values
inner_algo = pg.nlopt("slsqp"); inner_algo.maxeval = 5
algo = pg.algorithm(uda = pg.mbh(inner_algo, stop = 2, perturb = .2))
#algo.set_verbosity(10) # in this case this correspond to logs each 1 call to slsqp
pop = pg.population(prob = my_constrained_udp(), size = 1000 , seed=123)
pop.problem.c_tol = [1E-6] * 1 # get_nec + get_nic = 1, so multiplied by 1
pop = algo.evolve(pop)
res = pop.champion_x
return res
# running above optimization code for 1000 random initializations
for i in range(1000):
init_val = np.array([random.uniform(0, 1) for k in range(node_len)])
if perturb:
for k,v in perturb.items(): init_val[k] = v #setting perturbations in initial values
res = dyn_optimGMO(matL ,node_len ,init_val) # this function is defined here only
inits.append(init_val); results.append(res)
编辑 1：
正如@Ananda 在下面所建议的，我对 objective 函数进行了更改，将 运行ning 时间减少了近 7 倍。我已经使用 python multiprocessing 模块将代码重写为 运行 此代码超过 1000 iterations 。下面是我的新代码，我在其中尝试生成进程以并行处理迭代。由于我的系统只有 8 个线程，所以我将池大小限制为 5 个，因为 PyGMO 也使用内部并行化并且它也需要一些线程
import numpy as np
import pygmo as pg
matL = np.random.rand(300,300) ; node_len = 300
perturb = {12:1} # assign your perturb ID here
def optimizationFN(var):
results = []
inits = var[0]; perturb = var[1]
class my_constrained_udp:
def fitness(self, x):
obj1 = x[None,:] @ matL @ x[:,None] # @ is mat multiplication operator
ce1 = np.sum(inits) - np.sum(x)
return [obj1, ce1]
def n_objs(self): # no of objectives
return 1
def get_nec(self): #no of equality constraints
return 1
def get_nic(self): #no of in-equality constraints
return 0
def get_bounds(self): #lower and upper bounds: use this to perturb the node
lowerB = np.array([0]*node_len); upperB = np.array([1]*node_len)
if perturb:
for k,v in perturb.items():
lowerB[k] = v; upperB[k] = v
return (lowerB,upperB)
def gradient(self, x):
return pg.estimate_gradient_h(lambda x: self.fitness(x), x)
def dyn_optimGMO(matL, node_len ,inits):
'''
perturb should be a dict of node index and value as 0 or 1. Type(node_index) = int
'''
if perturb:
for k,v in perturb.items(): inits[k] = v #setting perturbations in initial values
inner_algo = pg.nlopt("slsqp"); inner_algo.maxeval = 5
algo = pg.algorithm(uda = pg.mbh(inner_algo, stop = 2, perturb = .2))
#algo.set_verbosity(10) # in this case this correspond to logs each 1 call to slsqp
pop = pg.population(prob = my_constrained_udp(), size = 100 , seed=123)
pop.problem.c_tol = [1E-6] * 1 # get_nec + get_nic = 1, so multiplied by 1
pop = algo.evolve(pop)
res = pop.champion_x
return res
if perturb:
for k,v in perturb.items(): inits[k] = v #setting perturbations in initial values
res = dyn_optimGMO(matL ,node_len ,inits) # this function is defined here only
results.append(res)
return results
import time
st = time.time()
#1000 random initialisations
initial_vals = []
for i in range(1000): initial_vals.append(np.array([random.uniform(0, 1) for k in range(node_len)]))
initial_vals = np.array(initial_vals)
inp_val = []
for i in range(len(initial_vals)): inp_val.append([initial_vals[i],perturb])
#running evaluations
#eqVal = optimizationFN(initial_vals,perturb=perturb)
from multiprocessing import Pool
myPool = Pool(8)
data = myPool.map(optimizationFN,inp_val)
myPool.close(); myPool.join()
print('Total Time: ',round(time.time()-st,4))
这会在 1.13 hours 中执行全部 1000 次迭代。
不过，有没有其他可能让我可以让它更快？
在尝试并行化等之前，尝试弄清楚性能瓶颈到底是什么并尝试解决它。
如果您使用 line profiler、
来分析您的适应度函数
Line # Hits Time Per Hit % Time Line Contents
==============================================================
28 @profile
29 def fitness(self, x):
30 96105 3577978.0 37.2 9.5 matA = np.matrix(x)
31 96105 8548353.0 88.9 22.7 obj1 = matA.dot(matL).dot(matA.T)[0,0] #this is int value
32 96105 25328835.0 263.6 67.4 ce1 = sum(init_val) - sum(x)
33 96105 121800.0 1.3 0.3 return [obj1, ce1]
如您所见，大部分时间花在 dot 和 sum 函数上，还有大量时间花在创建 matA.
上
我会像这样重写函数 -
def fitness(self, x):
obj1 = x[None, :] @ matL @ x[:, None]
ce1 = np.sum(init_val) - np.sum(x)
return [obj1, ce1]
如果你分析这个函数你可以看到，
Line # Hits Time Per Hit % Time Line Contents
==============================================================
20 @profile
21 def fitness(self, x):
22
23 77084 3151649.0 40.9 48.9 obj1 = x[None, :] @ matL @ x[:, None]
24 77084 3214012.0 41.7 49.9 ce1 = np.sum(init_val) - np.sum(x)
25
26 77084 79439.0 1.0 1.2 return [obj1, ce1]
完整功能每次点击的总时间从大约 380 下降到 80。
建议不再使用 np.matrix 方法并将弃用。并且使用原生 python sum 而不是 np.sum 会降低很多性能。
在我的机器上，它使性能从 33 sec/it 提高到 6 sec/it 次。性能提升约 5 倍。
Q : "Is there any way to parallelize this code and make it faster for any number of iterations?"
是的。如果尝试“从外部”numba.jit() 代码（由于 Numba 编译警告的原因而失败），您可以求助于分发关于上述 1k+ 独立初始化的批次的部分，并让计算这些并行并在之后收集结果。
这样做的好处是性能瞬间提高 1000 倍，并且可以进一步扩展。
如果您使用的是大约 1k+ 个节点的大学集群，您的 1k 批次计算可以在大约相同的时间内产生结果，而 solo-运行 将执行 1k 长序列中的第一个（此处的通信成本可以忽略不计，请参阅 ）。
我正在尝试使用 PyGMO package 求解非线性优化
优化 class 是单独定义的，然后通过单独的函数 dyn_optimGMO 调用。必须为变量（节点）定义的 1000 个随机初始向量完成并保存此优化 inits ( or init_val)
使用 timeit 模块我发现每次迭代大约需要 17 seconds 才能完成。这意味着 1000 iterations 大约需要 5 hours 。这是一个非常大的时间。
如果我必须对 20 perturb 个节点重复此操作，那么总迭代次数将达到 200000，这将花费上面计算的线性时间。
我尝试通过使用 python multiprocessing 模块为 20 个扰动节点中的每一个并行化每组 1000 次迭代来解决这个问题。但这于事无补。
我也尝试过使用 Numba jit 函数，但它们无法识别 pyGMO 模块，因此失败了。
有什么方法可以并行化此代码并使其在任意次数的迭代中都更快？
如果我的问题足够清楚，请告诉我，如果不够清楚，我会根据需要添加详细信息。
import numpy as np
import pygmo as pg
matL = np.random.rand(300,300) ; node_len = 300
inits = []; results = []
perturb = {25:0} #setting a random node, say, node 25 to 0
class my_constrained_udp:
def __init__(self):
pass
def fitness(self, x):
matA = np.matrix(x)
obj1 = matA.dot(matL).dot(matA.T)[0,0] #this is int value
ce1 = sum(init_val) - sum(x)
return [obj1, ce1]
def n_objs(self): # no of objectives
return 1
def get_nec(self): #no of equality constraints
return 1
def get_nic(self): #no of in-equality constraints
return 0
def get_bounds(self): #lower and upper bounds: use this to perturb the node
lowerB = np.array([0]*node_len); upperB = np.array([1]*node_len)
if perturb:
for k,v in perturb.items():
lowerB[k] = v; upperB[k] = v
return (lowerB,upperB)
def gradient(self, x):
return pg.estimate_gradient_h(lambda x: self.fitness(x), x)
def dyn_optimGMO(matL, node_len ,init):
if perturb:
for k,v in perturb.items(): init_val[k] = v #setting perturbations in initial values
inner_algo = pg.nlopt("slsqp"); inner_algo.maxeval = 5
algo = pg.algorithm(uda = pg.mbh(inner_algo, stop = 2, perturb = .2))
#algo.set_verbosity(10) # in this case this correspond to logs each 1 call to slsqp
pop = pg.population(prob = my_constrained_udp(), size = 1000 , seed=123)
pop.problem.c_tol = [1E-6] * 1 # get_nec + get_nic = 1, so multiplied by 1
pop = algo.evolve(pop)
res = pop.champion_x
return res
# running above optimization code for 1000 random initializations
for i in range(1000):
init_val = np.array([random.uniform(0, 1) for k in range(node_len)])
if perturb:
for k,v in perturb.items(): init_val[k] = v #setting perturbations in initial values
res = dyn_optimGMO(matL ,node_len ,init_val) # this function is defined here only
inits.append(init_val); results.append(res)
编辑 1：
正如@Ananda 在下面所建议的，我对 objective 函数进行了更改，将 运行ning 时间减少了近 7 倍。我已经使用 python multiprocessing 模块将代码重写为 运行 此代码超过 1000 iterations 。下面是我的新代码，我在其中尝试生成进程以并行处理迭代。由于我的系统只有 8 个线程，所以我将池大小限制为 5 个，因为 PyGMO 也使用内部并行化并且它也需要一些线程
import numpy as np
import pygmo as pg
matL = np.random.rand(300,300) ; node_len = 300
perturb = {12:1} # assign your perturb ID here
def optimizationFN(var):
results = []
inits = var[0]; perturb = var[1]
class my_constrained_udp:
def fitness(self, x):
obj1 = x[None,:] @ matL @ x[:,None] # @ is mat multiplication operator
ce1 = np.sum(inits) - np.sum(x)
return [obj1, ce1]
def n_objs(self): # no of objectives
return 1
def get_nec(self): #no of equality constraints
return 1
def get_nic(self): #no of in-equality constraints
return 0
def get_bounds(self): #lower and upper bounds: use this to perturb the node
lowerB = np.array([0]*node_len); upperB = np.array([1]*node_len)
if perturb:
for k,v in perturb.items():
lowerB[k] = v; upperB[k] = v
return (lowerB,upperB)
def gradient(self, x):
return pg.estimate_gradient_h(lambda x: self.fitness(x), x)
def dyn_optimGMO(matL, node_len ,inits):
'''
perturb should be a dict of node index and value as 0 or 1. Type(node_index) = int
'''
if perturb:
for k,v in perturb.items(): inits[k] = v #setting perturbations in initial values
inner_algo = pg.nlopt("slsqp"); inner_algo.maxeval = 5
algo = pg.algorithm(uda = pg.mbh(inner_algo, stop = 2, perturb = .2))
#algo.set_verbosity(10) # in this case this correspond to logs each 1 call to slsqp
pop = pg.population(prob = my_constrained_udp(), size = 100 , seed=123)
pop.problem.c_tol = [1E-6] * 1 # get_nec + get_nic = 1, so multiplied by 1
pop = algo.evolve(pop)
res = pop.champion_x
return res
if perturb:
for k,v in perturb.items(): inits[k] = v #setting perturbations in initial values
res = dyn_optimGMO(matL ,node_len ,inits) # this function is defined here only
results.append(res)
return results
import time
st = time.time()
#1000 random initialisations
initial_vals = []
for i in range(1000): initial_vals.append(np.array([random.uniform(0, 1) for k in range(node_len)]))
initial_vals = np.array(initial_vals)
inp_val = []
for i in range(len(initial_vals)): inp_val.append([initial_vals[i],perturb])
#running evaluations
#eqVal = optimizationFN(initial_vals,perturb=perturb)
from multiprocessing import Pool
myPool = Pool(8)
data = myPool.map(optimizationFN,inp_val)
myPool.close(); myPool.join()
print('Total Time: ',round(time.time()-st,4))
这会在 1.13 hours 中执行全部 1000 次迭代。
不过，有没有其他可能让我可以让它更快？
在尝试并行化等之前，尝试弄清楚性能瓶颈到底是什么并尝试解决它。
如果您使用 line profiler、
来分析您的适应度函数Line # Hits Time Per Hit % Time Line Contents
==============================================================
28 @profile
29 def fitness(self, x):
30 96105 3577978.0 37.2 9.5 matA = np.matrix(x)
31 96105 8548353.0 88.9 22.7 obj1 = matA.dot(matL).dot(matA.T)[0,0] #this is int value
32 96105 25328835.0 263.6 67.4 ce1 = sum(init_val) - sum(x)
33 96105 121800.0 1.3 0.3 return [obj1, ce1]
如您所见，大部分时间花在 dot 和 sum 函数上，还有大量时间花在创建 matA.
我会像这样重写函数 -
def fitness(self, x):
obj1 = x[None, :] @ matL @ x[:, None]
ce1 = np.sum(init_val) - np.sum(x)
return [obj1, ce1]
如果你分析这个函数你可以看到，
Line # Hits Time Per Hit % Time Line Contents
==============================================================
20 @profile
21 def fitness(self, x):
22
23 77084 3151649.0 40.9 48.9 obj1 = x[None, :] @ matL @ x[:, None]
24 77084 3214012.0 41.7 49.9 ce1 = np.sum(init_val) - np.sum(x)
25
26 77084 79439.0 1.0 1.2 return [obj1, ce1]
完整功能每次点击的总时间从大约 380 下降到 80。
建议不再使用np.matrix 方法并将弃用。并且使用原生 python sum 而不是 np.sum 会降低很多性能。
在我的机器上，它使性能从 33 sec/it 提高到 6 sec/it 次。性能提升约 5 倍。
Q : "Is there any way to parallelize this code and make it faster for any number of iterations?"
是的。如果尝试“从外部”numba.jit() 代码（由于 Numba 编译警告的原因而失败），您可以求助于分发关于上述 1k+ 独立初始化的批次的部分，并让计算这些并行并在之后收集结果。
这样做的好处是性能瞬间提高 1000 倍，并且可以进一步扩展。
如果您使用的是大约 1k+ 个节点的大学集群，您的 1k 批次计算可以在大约相同的时间内产生结果，而 solo-运行 将执行 1k 长序列中的第一个（此处的通信成本可以忽略不计，请参阅