python中一个三维优化问题的并行化

Question

我正在用scipy.optimize处理一个优化问题，该问题旨在计算一些3D地图。给定3个真实数据体(vol0, vol1, vol2)，我的目标是通过将体素强度的某些函数按体素方式拟合到其模型来估计3个参数(map_0, map_1, map_2)的贴图

到目前为止，我的想法是：

import scipy
import numpy as np
from scipy.optimize import minimize

def objective (x,arg0,arg1):

    vol_model0 = someFun( arg0[0], arg1, ... ) # *model value for arg0 which needs arg0[0] and arg1*
    vol_model1 = someFun( arg0[1], arg1, ... ) # *model value for arg0 which needs arg0[1] and arg1*
    vol_model2 = someFun( arg0[2], arg1, ... ) # *model value for arg0 which needs arg0[2] and arg1*

    RSS = np.sum( [ np.power( ( arg0[0] - vol_model0 ), 2 )
                  + np.power( ( arg0[1] - vol_model1 ), 2 )
                  + np.power( ( arg0[2] - vol_model2 ), 2 ) ]
                  )
    return RSS


arg1 = [1, 2, 3, 4]

vol0 = 5* np.zeros([100,100,100])
vol1 = 3* np.zeros([100,100,100])
vol2 = 4* np.zeros([100,100,100])

map_0 = np.zeros([100,100,100])
map_1 = np.zeros([100,100,100])
map_2 = np.zeros([100,100,100])

x0    = [5, 5, 5] 
bnds  = ( (1,10), (1, 10), (1, 10) )

for i0 in range(0,100):
    for i1 in range(0,100):
        for i2 in range(0,100):
            arg0 = [ vol0[i0,i1,i2], \
                     vol1[i0,i1,i2],  \
                     vol2[i0,i1,i2]    \
                     ]
            res = minimize(objective, x0, args = (arg0,arg1), bounds = bnds)
            map_0[i0,i1,i2], \
            map_1[i0,i1,i2],   \
            map_2[i0,i1,i2] = res.x

我的问题是：
考虑到这个约束优化问题，有没有办法让整个过程更快？嵌套的for循环需要很长时间才能完成。 我想知道是否有一种方法可以并行处理这个问题，或者使它更快

Answer 1

问题更多地与处理的效率有关，而不是与使用任何并发处理有关

在进行任何进一步的步骤之前，让我提出几个步骤：

1）
如果将完成1E6调用，请优化所有函数开销-因为保存的每个[us]将在终点线上保存另一个完整的[s]

2）
避免任何不相关的开销（为什么要让python创建一个变量，维护一个变量名，如果它永远不会被重复使用的话？你只需支付成本，而不会从中获得任何好处）

3）
尽最大努力提高性能和；如果可能的话，numba编译objective()函数的代码，
在这里，从{}运行时中删除的每一个{}都将进一步节省{}-在终点线上的重复次数的倍…

值得这样做，不是吗

---

即时奖励：
便宜的&；低垂果实：
矢量化&；充分利用Numpy工具的威力

>>> aClk.start(); _ = np.sum( # _______________________________ AVOID I/O on printing here by assignment into _ var
                              [ np.power( ( arg0[0] - v0 ), 2 ),
                                np.power( ( arg0[1] - v1 ), 2 ),
                                np.power( ( arg0[2] - v2 ), 2 )
                                ] #___________________________ LIST-constructor overheads
                              ); aClk.stop() # _______________.stop() the [us]-clock
225 [us]
164 [us]
188
157
175
199
201
201
171 [us]

>>> RSS = _
>>> RSS
0.16506628505715615

---

交叉验证方法的正确性：

>>> arg0_dif     = arg0.copy()
>>> vV           = np.ones( 3 ) * 1.23456789
>>> arg0_dif[:] -= vV[:]
>>> np.dot( arg0_dif, arg0_dif.T )
0.16506628505715615

---

+5x ~ +7x更快的代码：

>>> aClk.start(); _ = np.dot( arg0_dif, arg0_dif.T ); aClk.stop()
39
38
28
28
27 [us] ~ +5x ~ +7x FASTER just by using more efficient way to produce the RSS calculations

---

下一步：

"""
@numba.jit( signature =    'f8[:], f8[:], i8[:]', nogil = True, ... ) """
def aMoreEfficientObjectiveFUN( x,  arg0,  arg1 ):
    ###############################################################################################
    # BEST 2 VECTORISE THE someFun() into someFunVECTORISED(),
    #                      so as to work straight with arg0-vector and return vector of differences
    #rg0[:]   -= someFunVECTORISED( arg0,    arg1, ... ) ##########################################
    arg0[0]   -= someFun(           arg0[0], arg1, ... ) # *model value for arg0 which needs arg0[0] and arg1*
    arg0[1]   -= someFun(           arg0[1], arg1, ... ) # *model value for arg0 which needs arg0[0] and arg1*
    arg0[2]   -= someFun(           arg0[2], arg1, ... ) # *model value for arg0 which needs arg0[0] and arg1*

    return  np.dot( arg0, arg0.T ) # just this change gets ~ 5x ~ 7x FASTER processing ............

---

正确完成这一点后，剩下的就很简单了：

在一种情况下，修改后的overhead-strict Amdahl's Law的所有附加间接成本仍然可以证明支付以下所有附加成本是合理的：
+SER/DES（在发送参数之前进行酸洗）
+XFER-main2process
+SER/DES（在发送参数之前进行酸洗）
+XFER-process2main，
^{可以帮助您拆分1E6，完全独立的调用并将结果重新收集回map_tensor[:,i0,i1,i2]
~( map_0[i0,i1,i2], map_1[i0,i1,i2], map_2[i0,i1,i2] )

在一种情况下，所有附加开销成本都不合理，保持工作流程不在几个进程之间分配，如果确实好奇并热衷于实验，可以尝试在基于线程的并发处理中运行此过程，如无GIL处理（转义，在numpy-不需要持有GIL锁的部件，GIL驱动的重新{{}-ISSION）的成本可能会在重叠并发处理上表现出一些延迟掩蔽。体内测试将证明这一点，或者表明，在python-as-is生态系统中这样做没有额外的优势

---

结语：

那太好了&；公平地告诉我们您对运行时改进的反馈，部分来自代码效率（从最初的~ 27 minutes（原样）
在智能np.dot()技巧
+在完全向量化的objective( x, arg0, arg1 )之后，其中someFun()可以处理向量，智能在@numba.jit( ... )之后返回np.dot()装饰处理，如果这样的结果表明自己进一步更有效