Scipy 使用分组依据 pandas 数据框最小化

Scipy minimize with pandas dataframe with group by

我有一个数据框(下面的示例 df)并试图最小化它的成本函数。

GrpId = ['A','A','A','A','A','A','B','B','B','B','B','B','B']
col1 = [69.1,70.5,71.4,72.8,73.2,74.2,208.0,209.2,210.2,211.0,211.2,211.7,212.5]
col2 = [2,3.1,1.1,2.1,6.0,1.1,1.2,1.3,3.1,2.9,5.0,6.1,3.2]
d = {'GrpId':GrpId,'col1':col1,'col2':col2}

df1 = pd.DataFrame(d)

下面是最小化和成本函数。

col1_const=[0,0,0,0,60.0,0,0,0]
col2_const=[0,0,0,0,0,100.0,0,0]

def main(type1,type2,type3,df):
    vall0=[type1,type2,type3]
    res=minimize(cost_fun, vall0, args=(df), method = 'SLSQP', tol=0.01)

    [type1,type2,type3]=res.x

    return type1,type2,type3

def cost_fun(v, df):

    df['col1_res'][i] = np.where((df['col1'][i]!=np.nan), ((1/0.095)*(np.sqrt(df['col1'][i])-np.sqrt(col1_const[4]*(0.1*v[1]+v[2])**2)))**2 ,0)
    df['col2_res'][i] = np.where((df['col2'][i]!=np.nan), ((1/0.12)*(np.sqrt(df['col2'][i])-np.sqrt(col2_const[5]*(0.1*v[0]+v[2])**2)))**2 ,0)   
    
    res=0.5*np.sqrt(df['col1_res'][i]+df['col2_res'][i])

    return res

然后我在循环中迭代这个函数,如下所示,它工作但需要大量时间和内存,

df1['type1']=np.nan
df1['type2']=np.nan
df1['type3']=np.nan
df1['col3']=np.nan
df1['col1_res']=np.nan
df1['col2_res']=np.nan

for i in range(len(df1.GrpId)):
    if i==0:
        df1['type1'][i], df1['type2'][i], df1['type3'][i]= main(0.125, 0.125, 0.125,df1)
    else:
        df1['type1'][i], df1['type2'][i], df1['type3'][i]= main(df1['type1'][i-1], df1['type2'][i-1], df1['type3'][i-1],df1)
    df1['col3'][i]=df1['type1'][i]+df1['type2'][i]

请注意,我有更大的数据框,有更多的行和列,对于这个问题,我刚刚创建了一个示例 code/case。

我的问题是,

  1. 如何在不迭代的情况下做同样的事情
  2. col1_const[4] 值将根据组(按 GrpId 分组)更改 - 我有另一个函数来计算每个组的 col1_const[4] 值。在这种情况下,如何按组将此值传递给 cost_fun。

首先,我认为没有必要在 objective 函数中检查 != np.nan。相反,您可以清理数据框并将所有 np.nan 替换为零。 objective 函数在优化例程中被多次调用,因此应尽可能高效快速地编写。因此,我们删除了 np.where 的调用。另请注意,依赖索引变量 i 在外部范围内已知的事实是不好的做法,并且会使代码难以阅读。我会推荐这样的东西:

col1 = df1.col1.values[df1.col1.values != np.nan]
col2 = df1.col2.values[df1.col2.values != np.nan]
col1_const = np.array([0,0,0,0,60.0,0,0,0])
col2_const = np.array([0,0,0,0,0,100.0,0,0])

def cost_fun1(v, *args):
    i, col1, col2, col1_const, col2_const = args
    col1_res = ((1/0.095)*(np.sqrt(col1[i]) - np.sqrt(col1_const[4]*(0.1*v[1]+v[2])**2)))**2
    col2_res = ((1/0.12)*(np.sqrt(col2[i]) - np.sqrt(col2_const[5]*(0.1*v[0]+v[2])**2)))**2
    return 0.5*np.sqrt(col1_res + col2_res)

接下来,更重要的是,您要解决多个优化问题,而不是解决一个 large-scale 优化问题。从数学上讲,因为您的 objective 函数保证为正,所以您可以按照与 相同的方式重新表述您的问题。然后,cost_fun2 基本上 returns 所有索引 i 的所有 cost_fun1 的总和。使用一点重塑魔法,函数看起来几乎相同:

def cost_fun2(vv, *args):
    col1, col2, col1_const, col2_const = args
    v = vv.reshape(3, col1.size)
    col1_res = ((1/0.095)*(np.sqrt(col1) - np.sqrt(col1_const[4]*(0.1*v[1]+v[2])**2)))**2
    col2_res = ((1/0.12)*(np.sqrt(col2) - np.sqrt(col2_const[5]*(0.1*v[0]+v[2])**2)))**2
    return np.sum(0.5*np.sqrt(col1_res + col2_res))

然后,我们简单地解决问题,然后将解决方案的值写入数据框:

from scipy.optimize import minimize

# initial guess
x0 = np.ones(3*col1.size)

# solve the problem
res = minimize(lambda vv: cost_fun2(vv, col1, col2, col1_const, col2_const), x0=x0, method="trust-constr")

# write to dataframe
type1_vals, type2_vals, type3_vals = np.split(res.x, 3)
df1['type1'] = type1_vals
df1['type2'] = type2_vals
df1['type3'] = type3_vals

如果数据框中需要 col1_rescol2_res,可以直接修改 objective 函数。

最后但同样重要的是,根据数据帧的大小,强烈建议将精确的 objective 梯度传递给 scipy.optimize.minimize 以获得良好的收敛性能。目前,梯度是通过有限差分来近似的,这很慢并且容易出现舍入误差。