python中3D优化问题的并行化

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() 函数的代码，
这里每个 [us] 从 someFun( aFloat64, anInt64[], ... ) 运行时中删除将进一步节省您 ~ 3 [s]-times-the-number-of-minimize()-er-在终点线上重新迭代...

值得这样做，不是吗？

---

即时奖励：
廉价且容易实现的果实：
矢量化和利用 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，
multiprocessing 可以帮助您拆分 1E6，完全独立的调用并将结果重新收集到map_tensor[:,i0,i1,i2]
~ ( map_0[i0,i1,i2], map_1[i0,i1,i2], map_2[i0,i1,i2] )

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

---

结语：

如果您向我们提供有关运行时改进的反馈，部分来自代码效率（从最初的 ~ 27 minutes（原样）
-- 在智能 @ 之后987654342@ 技巧
-- +在完全矢量化的objective( x, arg0, arg1 ) 之后，其中someFun() 可以处理矢量并智能返回 np.dot() 本身
-- +在 @numba.jit( ... ) 修饰处理之后，如果这样的结果表明自己更加高效