Графики Scipy optimize. Минимизировать результаты сходимости на каждой итерации? - PullRequest
Новая
       21

Графики Scipy optimize. Минимизировать результаты сходимости на каждой итерации?

0 голосов
/ 07 октября 2018

Я хотел бы выполнить некоторые тесты для моей процедуры оптимизации, используя scipy.optimize.minimize, в частности, построить график сходимости (или, скорее, объективной функции) каждой итерации по нескольким тестам.

Предположим, у меня есть следующая линейно ограниченная задача квадратичной оптимизации:

свести к минимуму: x_i Q_ij x_j + a | x_i |

подчиняется: sum (x_i) =1

Я могу кодировать это как: 

def _fun(x, Q, a):
    c = np.einsum('i,ij,j->', x, Q, x)
    p = np.sum(a * np.abs(x))
    return c + p
def _constr(x):
    return np.sum(x) - 1

И я буду реализовывать оптимизацию в scipy как: 

x_0 = # some initial vector
x_soln = scipy.optimise.minimize(_fun, x_0, args=(Q, a), method='SLSQP', 
                                 constraints={'type': 'eq', 'fun': _constr})

Я вижу, что тамcallback аргумент , но принимает только один аргумент значений параметров на каждой итерации .Как я могу использовать это в более эзотерическом случае, когда у меня могут быть другие аргументы, которые должны быть предоставлены моей функции обратного вызова?

1 Ответ

0 голосов
/ 07 октября 2018

Я решил это, используя общий объект кэша обратного вызова, на который ссылались каждый раз из моей функции обратного вызова.Допустим, вы хотите сделать 20 тестов и построить целевую функцию после каждой итерации на одном и том же графике.Вам понадобится внешний цикл для запуска 20 тестов, но мы создадим его позже.

Сначала давайте создадим класс, который будет хранить для нас все значения целевой функции итерации, а также пару дополнительных кусочков: 

class OpObj(object):
    def __init__(self, Q, a):
        self.Q, self.a = Q, a
        rv = np.random.rand()
        self.x_0 = np.array([rv, (1-rv)/2, (1-rv)/2])
        self.f = np.full(shape=(500,), fill_value=np.NaN)
        self.count = 0
    def _fun(self, x):
        return _fun(x, self.Q, self.a)

Также позволяет добавить функцию обратного вызова, которая управляет этим классом obj.Не беспокойтесь, у него есть более одного аргумента, так как мы исправим это позже.Просто убедитесь, что первый параметр - это переменные решения. 

def cb(xk, obj=None):
    obj.f[obj.count] = obj._fun(xk)
    obj.count += 1

Все, что он делает, это использует функции и значения объекта для обновления, подсчитывая количество итераций каждый раз.Эта функция будет вызываться после каждой итерации.

Собирая все это вместе, все, что нам нужно, это еще две вещи: 1) некоторая математическая привязка для построения графика и исправление обратного вызова, чтобы иметь только один аргумент.Мы можем сделать это с помощью декоратора, который является именно тем, что делает частично функция functools.Возвращает функцию с меньшим количеством аргументов, чем у оригинала.Итак, окончательный код выглядит так: 

import matplotlib.pyplot as plt
import scipy.optimize as op
import numpy as np
from functools import partial

Q = np.array([[1.0, 0.75, 0.45], [0.75, 1.0, 0.60], [0.45, 0.60, 1.0]])
a = 1.0

def _fun(x, Q, a):
    c = np.einsum('i,ij,j->', x, Q, x)
    p = np.sum(a * np.abs(x))
    return c + p
def _constr(x):
    return np.sum(x) - 1

class OpObj(object):
    def __init__(self, Q, a):
        self.Q, self.a = Q, a
        rv = np.random.rand()
        self.x_0 = np.array([rv, (1-rv)/2, (1-rv)/2])
        self.f = np.full(shape=(500,), fill_value=np.NaN)
        self.count = 0
    def _fun(self, x):
        return _fun(x, self.Q, self.a)

def cb(xk, obj=None):
    obj.f[obj.count] = obj._fun(xk)
    obj.count += 1

fig, ax = plt.subplots(1,1)
x = np.linspace(1,500,500)
for test in range(20):
    op_obj = OpObj(Q, a)
    x_soln = op.minimize(_fun, op_obj.x_0, args=(Q, a), method='SLSQP',
                         constraints={'type': 'eq', 'fun': _constr},
                         callback=partial(cb, obj=op_obj))
    ax.plot(x, op_obj.f)

ax.set_ylim((1.71,1.76))
plt.show()

enter image description here

...