Scipy.optimize.minimize(method=‘slsqp’) Memory Issues When Outside The Bounds

I am working with scipy.optimize.minimize(method=‘SLSQP’), the function and constraints are interpolated with scipy.interpolate.LinearNDInterpolator. The start values are rand

Solution 1:

My own solution up to now is not to use the scipy.interpolate.LinearNDInterpolator but to replace it with several one-dimensional interpolations. For a 3D Problem and a cubic interpolation with scipy.interpolate.UnivariateSpline I need 64 data points on a rectangular grid. In the first step I only interpolate (16 times) in coordinate direction 1 and reduce my problem to a 2D-one with 16 points. In the next step I interpolate (4 times) in coordinate direction 2 and reduce further to a 1D promlem with 4 points left. The last step is one 1D interpolation in coordinate direction 3. With this procedure scipy.optimize.minimize(method=‘SLSQP’) does not run into memory issues even when outside the bounds. The interpolation code is the following:

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

defsort2lists(list1=[],list2=[]):
    '''
    Lists must be equal size
    sorts after list1
    '''
    list1list2 = [(list1[i1],list2[i1]) for i1 inrange(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


defsortdictkeysandreturnvalues(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

defmeshgrid2(*arrs):
    iflen(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 = 1for s in lens:
    sz*=s

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


defarr2points(*arrs):
    '''
    arrs = ([list1], [list2], ...)
    ans = [[x1,y1,z1...],[x1,y1,z2,...],...,[x1,y2,z1,...],...[x2,y1,z1]
    '''iflen(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 = 1for s in lens:
    sz*=s

    ans = [[] for i1 inrange(sz)]    

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

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

#defN1Dsplines(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 inrange(k+1):
        points = [i2 for i2 in pointsvaluesdict.keys() if i2[i1] notin 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 inrange(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 inrange(len(pointsgrid))}
    return griddict
#

Test 4D quadratic

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)

Solution 2:

Still, why was your interpolation going outside your defined bounds for your variables in the first place?

It seems like you defined your max and min bounds for your variables.

I know that depending on the version of scipy I use, sometimes the minimize function will not accept my bounds, unless they are given in floats ( it's a really terrible bug ).

You could try redefining your bounds as

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

This usually works for me.

Solution 3:

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

solved the problem completely