scipy.optimize.minimize（method ='SLSQP'）超出界限时的内存问题

Question

我正在使用scipy.optimize.minimize(method=‘SLSQP’)，函数和约束使用scipy.interpolate.LinearNDInterpolator进行插值。起始值是边界内的随机数。

我一起工作：

1. SciPy的0.13.3
2. 用Cython 0.20.1

的优化运行，有时并给出一个合理的结果，但有时优化开始要求巨额的内存高达20GB，然后我的电脑停止工作。这总是出现在边界之外的值上。

难道scipy.interpolate.LinearNDInterpolator不能与scipy.optimize.minimize(method=‘SLSQP’)一起使用吗？在边界之外，我没有模拟数据，所以插值给出了fill_Value = 0或fill_value = 1e10。当我与scipy.optimize.fmin_slsqp工作发生

相同的行为

不幸的是我的代码是非常大的，但与这组数据我总是得到的内存问题：

######################################### 
###Memory Leak scipy.optimize.minimize### 
######################################### 
import numpy as np 
from scipy.optimize import minimize 
from scipy.interpolate import LinearNDInterpolator 
def objfun(x): 
    print x 
    return x[1] 

points = np.array([[ 0.00000000e+00, 0.00000000e+00],[ 0.00000000e+00, 1.00334000e+00],[ 0.00000000e+00, 2.00669000e+00],[ 7.07700000e+02, 0.00000000e+00],[ 7.07700000e+02, 1.00334000e+00],[ 7.07700000e+02, 2.00669000e+00],[ 1.56890000e+03, 0.00000000e+00],[ 1.56890000e+03, 1.00334000e+00],[ 1.56890000e+03, 2.00669000e+00],[ 2.50080000e+03, 0.00000000e+00],[ 2.50080000e+03, 1.00334000e+00],[ 2.50080000e+03, 2.00669000e+00],[ 3.47090000e+03, 0.00000000e+00],[ 3.47090000e+03, 1.00334000e+00],[ 3.47090000e+03, 2.00669000e+00],[ 4.46380000e+03, 0.00000000e+00],[ 4.46380000e+03, 1.00334000e+00],[ 4.46380000e+03, 2.00669000e+00],[ 5.47130000e+03, 0.00000000e+00],[ 5.47130000e+03, 1.00334000e+00],[ 5.47130000e+03, 2.00669000e+00],[ 6.48890000e+03, 0.00000000e+00],[ 6.48890000e+03, 1.00334000e+00],[ 6.48890000e+03, 2.00669000e+00],[ 7.51360000e+03, 0.00000000e+00],[ 7.51360000e+03, 1.00334000e+00],[ 7.51360000e+03, 2.00669000e+00],[ 8.54350000e+03, 0.00000000e+00],[ 8.54350000e+03, 1.00334000e+00],[ 8.54350000e+03, 2.00669000e+00],[ 9.57740000e+03, 0.00000000e+00],[ 9.57740000e+03, 1.00334000e+00],[ 9.57740000e+03, 2.00669000e+00],[ 1.06143000e+04, 0.00000000e+00],[ 1.06143000e+04, 1.00334000e+00],[ 1.06143000e+04, 2.00669000e+00],[ 1.16535000e+04, 0.00000000e+00],[ 1.16535000e+04, 1.00334000e+00],[ 1.16535000e+04, 2.00669000e+00],[ 1.26945000e+04, 0.00000000e+00],[ 1.26945000e+04, 1.00334000e+00],[ 1.26945000e+04, 2.00669000e+00],[ 1.37369000e+04, 0.00000000e+00],[ 1.37369000e+04, 1.00334000e+00],[ 1.37369000e+04, 2.00669000e+00],[ 1.47804000e+04, 0.00000000e+00],[ 1.47804000e+04, 1.00334000e+00],[ 1.47804000e+04, 2.00669000e+00],[ 1.58248000e+04, 0.00000000e+00],[ 1.58248000e+04, 1.00334000e+00],[ 1.58248000e+04, 2.00669000e+00],[ 1.68698000e+04, 0.00000000e+00],[ 1.68698000e+04, 1.00334000e+00],[ 1.68698000e+04, 2.00669000e+00],[ 1.79153000e+04, 0.00000000e+00],[ 1.79153000e+04, 1.00334000e+00],[ 1.79153000e+04, 2.00669000e+00],[ 1.89612000e+04, 0.00000000e+00],[ 1.89612000e+04, 1.00334000e+00],[ 1.89612000e+04, 2.00669000e+00],[ 2.00074000e+04, 0.00000000e+00],[ 2.00074000e+04, 1.00334000e+00],[ 2.00074000e+04, 2.00669000e+00]]) 
values = np.array([ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,4.29730000e+01, 5.72947500e-01, -5.35464000e-01,9.11676000e+01, 1.31063500e+00, -1.05937500e+00,1.38660750e+02, 2.11484000e+00, -1.50850500e+00,1.84497000e+02, 2.96052000e+00, -1.88466000e+00,2.28622000e+02, 3.83846000e+00, -2.19702000e+00,2.71163000e+02, 4.74426500e+00, -2.45397000e+00,3.12274500e+02, 5.67547500e+00, -2.66222500e+00,3.52102000e+02, 6.63058000e+00, -2.82711000e+00,3.90774000e+02, 7.60858000e+00, -2.95286000e+00,4.28399500e+02, 8.60879000e+00, -3.04289000e+00,4.65074500e+02, 9.63071000e+00, -3.10001500e+00,5.00881500e+02, 1.06739850e+01, -3.12655500e+00,5.35893000e+02, 1.17383500e+01, -3.12444000e+00,5.70166500e+02, 1.28235000e+01, -3.09540500e+00,6.03760000e+02, 1.39293500e+01, -3.04082500e+00,6.36721500e+02, 1.50557500e+01, -2.96194500e+00,6.69093500e+02, 1.62026000e+01, -2.85982000e+00,7.00915000e+02, 1.73698000e+01, -2.73539500e+00,7.32222000e+02, 1.85573500e+01, -2.58950000e+00,7.63042500e+02, 1.97651000e+01, -2.42286000e+00]) 

S22_Adh1Ad_inter = LinearNDInterpolator(points,values,1e10) 
def Fsigcon(x): 
    rf1_int = x[1] 
    rf_eval=[] 
    x_eval=[] 
    interval = np.linspace(0,x[0],x[0]/0.01) 
    if interval.size == 0: 
     interval=np.array([x[0]]) 
    for xcoord in interval: 
     rf_eval.append(rf1_int) 
     x_eval.append(xcoord) 
    val_interp = S22_Adh1Ad_inter(rf_eval,x_eval) #'nearest' #'linear' #'cubic' 
    out = (val_interp.min()-39.45) 
    return out 

points = np.array([[ 0.00000000e+00, 0.00000000e+00],[ 0.00000000e+00, 1.99997000e-01],[ 0.00000000e+00, 4.00002000e-01],[ 7.07700000e+02, 1.39999000e-01],[ 7.07700000e+02, 3.39996000e-01],[ 1.56890000e+03, 8.00020000e-02],[ 1.56890000e+03, 2.79999000e-01],[ 2.50080000e+03, 1.99970000e-02],[ 2.50080000e+03, 2.20001000e-01],[ 2.50080000e+03, 1.90000200e+00],[ 3.47090000e+03, 1.60004000e-01],[ 3.47090000e+03, 3.60001000e-01],[ 4.46380000e+03, 9.99980000e-02],[ 4.46380000e+03, 3.00003000e-01],[ 5.47130000e+03, 4.00010000e-02],[ 5.47130000e+03, 2.39998000e-01],[ 5.47130000e+03, 3.00000000e+00],[ 6.48890000e+03, 1.80000000e-01],[ 6.48890000e+03, 3.79997000e-01],[ 7.51360000e+03, 1.20003000e-01],[ 7.51360000e+03, 3.20000000e-01],[ 8.54350000e+03, 5.99980000e-02],[ 8.54350000e+03, 2.60002000e-01],[ 9.57740000e+03, 0.00000000e+00],[ 9.57740000e+03, 1.99997000e-01],[ 9.57740000e+03, 4.00002000e-01],[ 1.06143000e+04, 1.39999000e-01],[ 1.06143000e+04, 3.39996000e-01],[ 1.16535000e+04, 8.00020000e-02],[ 1.16535000e+04, 2.79999000e-01],[ 1.26945000e+04, 1.99970000e-02],[ 1.26945000e+04, 2.20001000e-01],[ 1.26945000e+04, 1.90000200e+00],[ 1.37369000e+04, 1.60004000e-01],[ 1.37369000e+04, 3.60001000e-01],[ 1.47804000e+04, 9.99980000e-02],[ 1.47804000e+04, 3.00003000e-01],[ 1.58248000e+04, 4.00010000e-02],[ 1.58248000e+04, 2.39998000e-01],[ 1.58248000e+04, 3.00000000e+00],[ 1.68698000e+04, 1.80000000e-01],[ 1.68698000e+04, 3.79997000e-01],[ 1.79153000e+04, 1.20003000e-01],[ 1.79153000e+04, 3.20000000e-01],[ 1.89612000e+04, 5.99980000e-02],[ 1.89612000e+04, 2.60002000e-01],[ 2.00074000e+04, 0.00000000e+00],[ 2.00074000e+04, 1.99997000e-01],[ 2.00074000e+04, 4.00002000e-01]]) 

values = np.array([ 0.  , 0.  , 0.  , 0.010168, 0.010055, 0.046252,0.045731, 0.092687, 0.107056, 0.11196 , 0.19232 , 0.190859,0.29924 , 0.295611, 0.401297, 0.42018 , 0.450553, 0.564416,0.561699, 0.727387, 0.719631, 0.883825, 0.894486, 0.  ,1.087256, 1.084631, 1.298136, 1.287209, 1.507127, 1.505308,1.424393, 1.740491, 1.839568, 1.993769, 1.981605, 2.251336,2.238475, 2.330676, 2.511822, 2.723058, 2.803453, 2.792818,3.104855, 3.08533 , 3.29549 , 3.393902, 0.  , 3.721085,3.714504]) 

G_Adh1Ad_inter = LinearNDInterpolator(points,values,0) 
def Gcon(x): 
    val_interp = G_Adh1Ad_inter(x[1],x[0]) 
    out = (val_interp.min()-0.33) 
    return out 

cons = (
    {'type': 'ineq', 
    'fun' : Fsigcon}, 
    {'type': 'ineq', 
    'fun' : Gcon} 
    ) 
amin = 0.0 
amax = 3.0 
bounds=[(amin,amax),(0.0, 20007.400000000001)] 
a_start= 1.5343936873636999 
rf1_start= 6824.9659188661817 
res_int = minimize(objfun, [a_start,rf1_start],method='SLSQP',jac=None,bounds=bounds,constraints=cons,tol =1e-4,options={'iprint':2, 'disp': True , 'maxiter':1e2}) 

Answer 1

自己的解决方案达现在不是使用scipy.interpolate.LinearNDInterpolator，而是用几个一维插值替换它。 对于3D问题和使用scipy.interpolate.UnivariateSpline的三次插值我需要矩形网格上的64个数据点。在第一步中，我只在坐标方向1上插入（16次），并将我的问题减少到16点的2D-1。在接下来的步骤中，我在坐标方向2上插值（4次），并进一步减少到4D点。最后一步是在坐标方向3上进行一次1D插值。 通过此过程，即使在边界之外，scipy.optimize.minimize（method ='SLSQP'）也不会遇到内存问题。 插值代码如下：

#!/usr/bin/python 
# coding: utf8 
# Filename: N1Dsplines.py 
import numpy as np 
from scipy.interpolate import UnivariateSpline 
from scipy.interpolate import LinearNDInterpolator 

def sort2lists(list1=[],list2=[]): 
    ''' 
    Lists must be equal size 
    sorts after list1 
    ''' 
    list1list2 = [(list1[i1],list2[i1]) for i1 in range(len(list1))] 
    list1list2 = sorted(list1list2, key=lambda x:x[0]) 
    list1 = [i1[0] for i1 in list1list2] 
    list2 = [i1[1] for i1 in list1list2] 
    return list1,list2 


def sortdictkeysandreturnvalues(dictionary): 
    keyvalues = [[i1,dictionary[i1]] for i1 in dictionary] 
    keyvalues = sorted(keyvalues, key= lambda x: x[0]) 
    points = [i1[1] for i1 in keyvalues] 
    return points 

def meshgrid2(*arrs): 
    if len(arrs)==1: 
    arrs = tuple(arrs[0]) 

    # list with entries= length of the input arrays 
    lens = map(len, arrs) 
    # dim = number of arrays 
    dim = len(arrs) 

    # sz is the multiplication of the length of the *arrs 
    sz = 1 
    for s in lens: 
    sz*=s 

    ans = []  
    for i, arr in enumerate(arrs): 
    slc = [1]*dim 
    slc[i] = lens[i] 
    arr2 = np.asarray(arr).reshape(slc) 
    for j, sz in enumerate(lens): 
     if j!=i: 
     arr2 = arr2.repeat(sz, axis=j) 
    ans.append(arr2) 
    return tuple(ans) 


def arr2points(*arrs): 
    ''' 
    arrs = ([list1], [list2], ...) 
    ans = [[x1,y1,z1...],[x1,y1,z2,...],...,[x1,y2,z1,...],...[x2,y1,z1] 
    ''' 
    if len(arrs)==1: 
    arrs = tuple(arrs[0]) 
    # list with entries= length of the input arrays 
    lens = map(len, arrs) 
    # dim = number of arrays 
    dim = len(arrs) 

    # sz is the multiplication of the length of the *arrs 
    sz = 1 
    for s in lens: 
    sz*=s 

    ans = [[] for i1 in range(sz)]  

    for i1, arr in enumerate(arrs): 
    count = 0 
    sz2 = 1 
    for s in lens[i1+1:]: 
     sz2*=s 

    for i4 in range(sz/(lens[i1]*sz2)): 
     for i2 in arr: 
     for i3 in range(sz2): 
      ans[count].append(i2) 
      count+=1 
    return ans 

# 





def N1Dsplines(pointsvaluesdict={}, point =(), scheme = (1,0), k=3): 
    ''' 
    Geht nur für reguläre Grids nicht für scattered data points 
    pointsvaluesdict = {(x,y,z):value , (x,y,z):value ....} 
    point = (x,y,z) 
    k = interpolation order 
    scheme = (2,1,0) bedeutet erst wird in Richtung z=2 dann in Richtung y=1 und als letztes in Richtung x=0 interpoliert 
    ''' 

    #Find the closest Datapoints to the point coordinates 
    closest = {i1:[] for i1 in scheme} 
    for i1 in scheme:  
    for count in range(k+1): 
     points = [i2 for i2 in pointsvaluesdict.keys() if i2[i1] not in closest[i1]] 
     closest[i1].append(min(points, key=lambda x:abs(x[i1]-point[i1]))[i1]) 

    #Make pointgrid from closest points 
    liste = sortdictkeysandreturnvalues(closest) 
    pointsgrid = arr2points(liste) 
    pointsgrid = [tuple(i1) for i1 in pointsgrid] 
    valuesgrid = [pointsvaluesdict[i1] for i1 in pointsgrid] 
    griddict = {pointsgrid[i1]:valuesgrid[i1] for i1 in range(len(pointsgrid))} 

    #Get 1D Splines by letting scheme[>0]=const and vary scheme[0] 
    for i1 in scheme: 
    pointsinter = [] 
    valuesinter = [] 
    temp = closest.copy() 
    del temp[i1] 
    points = sortdictkeysandreturnvalues(temp) 
    points = arr2points(points) 
    for coords in points: 

     points2inter=[] 
     for i2 in pointsgrid: 
     i2red = list(i2) 
     del i2red[i1] 
     if coords == i2red: 
      points2inter.append(i2) 
     values2inter = [griddict[i2] for i2 in points2inter] 
     points2inter = [i2[i1] for i2 in points2inter] 
     points2inter,values2inter = sort2lists(list1=points2inter,list2=values2inter) 

     coords.insert(i1,point[i1]) 
     pointsinter.append(coords) 
     s = UnivariateSpline(points2inter, values2inter,k=k, s=0) 
     value = s(point[i1]) 
     valuesinter.append(value) 

     closest[i1]=[point[i1]] 

    pointsgrid = pointsinter 
    pointsgrid = [tuple(i1) for i1 in pointsgrid] 
    valuesgrid = valuesinter 
    griddict = {pointsgrid[i1]:valuesgrid[i1] for i1 in range(len(pointsgrid))} 
    return griddict 
# 

测试4D二次

x = np.array([1, 2, 3, 4, 5, 6, 7 ,8 ,9, 10]) 
y = x 
z= x 
w = x 

xx, yy, zz, ww = meshgrid2(x, y,z,w) 
v= xx**2 + yy**2 +zz**2 +ww**2 
xx = xx.flatten() 
yy = yy.flatten() 
zz = zz.flatten() 
ww = ww.flatten() 
v = v.flatten() 
pointsvalues = [((xx[i1],yy[i1],zz[i1], ww[i1]),v[i1]) for i1 in range(len(v))] 
pointsvaluesdict = {key: value for (key, value) in pointsvalues} 
pointsvaluesdict.keys() 
pointsvaluesdict.values() 

point = (5.5,4.5,3.5,2.5) 
5.5**2 + 4.5**2 + 3.5**2 +2.5**2== 69 
N1Dsplines(pointsvaluesdict=pointsvaluesdict, point =(5.5,4.5,3.5,2.5), scheme = (0,1,2,3), k=3) 

#Outside Bounds 
N1Dsplines(pointsvaluesdict=pointsvaluesdict, point =(0,0,0,0), scheme = (0,1,2,3), k=3) 

Answer 2

尽管如此，为什么你要插你定义的边界之外为您的变量摆在首位？

好像你定义你的max和min界的变量。

我知道，根据我使用的scipy的版本，有时候最小化函数不会接受我的界限，除非它们在浮点数中给出（这是一个非常糟糕的错误）。

你可以尝试重新定义你的边界作为

bnds = np.array([ [min_var1, max_var1], [min_var2, max_var2] ], dtype = float) 

这通常为我工作。

Answer 3

sudo -H pip install --upgrade scipy 
sudo -H pip install --upgrade numpy 

解决了这个问题完全