Scipy: Minimize viole les bornes données

Question

Je cherche à utiliser scipy.optimize.minimize avec des bornes simples a <= x <= b. Cependant, il arrive souvent que ma fonction cible soit évaluée juste à l'extérieur des bornes. À ma connaissance, cela se produit lorsque minimize tente de déterminer le gradient de la fonction cible à la frontière.

---

Exemple minimal:

import math
import numpy as np
from scipy.optimize import Bounds, minimize

constraint = Bounds([-1, -1], [1, 1], True)
def fun(x):
    print(x)
    return -math.exp(-np.dot(x,x))
result = minimize(fun, [-1, -1], bounds=constraint)

La sortie montre que le minimiseur saute au point [1,1] puis tente d'évaluer à [1,00000001, 1]:

[-1. -1.]
[-0.99999999 -1.        ]
[-1.         -0.99999999]
[-0.72932943 -0.72932943]
[-0.72932942 -0.72932943]
[-0.72932943 -0.72932942]
[-0.22590689 -0.22590689]
[-0.22590688 -0.22590689]
[-0.22590689 -0.22590688]
[1. 1.]
[1.00000001 1.        ]
[1.         1.00000001]
[-0.03437328 -0.03437328]
...

Évidemment, il n'y a aucun problème dans cet exemple, car fun peut être évalué aussi là. Mais ce n'est pas toujours le cas...

---

Dans mon problème réel, le minimum ne peut pas être sur la frontière et j'ai la solution de contournement facile d'ajouter un epsilon aux bornes. Mais on pourrait s'attendre à ce qu'il y ait une solution simple à ce problème qui fonctionne également si le minimum peut être à une frontière?

PS : Il serait étrange que je sois le premier à rencontrer ce problème -- désolé si cette question a été posée auparavant quelque part, mais je ne l'ai pas trouvée.

Answer 1

Comme discuté ici (merci @"Welcome to Stack Overflow" pour le commentaire qui m'y a dirigé), le problème est effectivement que la routine de gradient ne respecte pas les limites. J'en ai écrit une nouvelle qui fait le travail :

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

def gradient_respecting_bounds(bounds, fun, eps=1e-8):
    """bounds: liste de tuples (inférieur, supérieur)"""
    def gradient(x):
        fx = fun(x)
        grad = np.zeros(len(x))
        for k in range(len(x)):
            d = np.zeros(len(x))
            d[k] = eps if x[k] + eps <= bounds[k][1] else -eps
            grad[k] = (fun(x + d) - fx) / d[k]
        return grad
    return gradient

bounds = ((-1, 1), (-1, 1))
def fun(x):
    print(x)
    return -math.exp(-np.dot(x,x))
result = minimize(fun, [-1, -1], bounds=bounds,
                  jac=gradient_respecting_bounds(bounds, fun))

---

Notez que cela peut être un peu moins efficace, car fun(x) est maintenant évalué deux fois à chaque point. Cela semble inévitable, extrait pertinent de _minimize_lbfgsb dans lbfgsb.py :

if jac is None:
    def func_and_grad(x):
        f = fun(x, *args)
        g = _approx_fprime_helper(x, fun, epsilon, args=args, f0=f)
        return f, g
else:
    def func_and_grad(x):
        f = fun(x, *args)
        g = jac(x, *args)
        return f, g

Comme vous pouvez le voir, la valeur de f ne peut être réutilisée que par la fonction interne _approx_fprime_helper.