Стохастическая оптимизация в Python

Я пытаюсь объединить cvxopt (решатель оптимизации) и PyMC (сэмплер) для решения выпуклых задач стохастической оптимизации.

Для справки, установка обоих пакетов с помощью pip проста:

pip install cvxopt
pip install pymc

Оба пакета отлично работают независимо друг от друга. Вот пример того, как решить проблему LP с cvxopt:

# Testing that cvxopt works
from cvxopt import matrix, solvers

# Example from http://cvxopt.org/userguide/coneprog.html#linear-programming

c = matrix([-4., -5.])
G = matrix([[2., 1., -1., 0.], [1., 2., 0., -1.]])
h = matrix([3., 3., 0., 0.])
sol = solvers.lp(c, G, h)
# The solution sol['x'] is correct: (1,1)

Однако, когда я пытаюсь использовать его с PyMC (например, поставив распределение на один из коэффициентов), PyMC выдает ошибку:

import pymc as pm
import cvxopt

c1 = pm.Normal('c1', mu=-4, tau=.5**-2)

@pm.deterministic
def my_lp_solver(c1=c1):
    c = matrix([c1, -5.])
    G = matrix([[2., 1., -1., 0.], [1., 2., 0., -1.]])
    h = matrix([3., 3., 0., 0.])
    sol = solvers.lp(c, G, h)
    solution = np.array(sol['x'],dtype=float).flatten()
    return solution

m = pm.MCMC(dict(c1=c1, x=x))
m.sample(20000, 10000, 10)

Я получаю следующую ошибку PyMC:

<ipython-input-21-5ce2909be733> in x(c1)
     14 @pm.deterministic
     15 def x(c1=c1):
---> 16     c = matrix([c1, -5.])
     17     G = matrix([[2., 1., -1., 0.], [1., 2., 0., -1.]])
     18     h = matrix([3., 3., 0., 0.])

TypeError: invalid type in list

Почему? Есть ли способ заставить cvxopt красиво играть с PyMC?

Фон:

Если кому-то интересно, PyMC позволяет вам использовать любую функцию по вашему выбору. В этом конкретном случае функция, из которой мы делаем выборку, — это функция, которая отображает проблему LP в решение. Мы делаем выборку из этой функции, потому что наша задача ЛП содержит стохастические коэффициенты, поэтому нельзя просто применить готовый решатель ЛП.

В частности, в этом случае один выходной образец PyMC — это просто решение проблемы LP. Поскольку параметры задачи LP изменяются (в соответствии с выбранными вами распределениями), выходные выборки из PyMC будут другими, и мы надеемся получить апостериорное распределение.

Приведенное выше решение основано на этом ответе, единственная разница в том, что я надеюсь использовать настоящий общий решатель ( в этом случае cvxopt)


python scipy mathematical-optimization cvxopt pymc
person Amelio Vazquez-Reina    schedule 01.08.2014    source источник

Ответы (1)


arrow_upward
5
arrow_downward

Тип c1, созданный с помощью pm.Normal, — это numpy array, вам просто нужно удалить его и преобразовать в float(c1), после чего он отлично работает:

>>> @pm.deterministic
... def my_lp_solver(c1=c1):
...     c = matrix([float(c1), -5.])
...     G = matrix([[2., 1., -1., 0.], [1., 2., 0., -1.]])
...     h = matrix([3., 3., 0., 0.])
...     sol = solvers.lp(c, G, h)
...     solution = np.array(sol['x'],dtype=float).flatten()
...     return solution
... 
     pcost       dcost       gap    pres   dres   k/t
 0: -8.1223e+00 -1.8293e+01  4e+00  0e+00  7e-01  1e+00
 1: -8.8301e+00 -9.4605e+00  2e-01  1e-16  4e-02  3e-02
 2: -9.0229e+00 -9.0297e+00  2e-03  2e-16  5e-04  4e-04
 3: -9.0248e+00 -9.0248e+00  2e-05  3e-16  5e-06  4e-06
 4: -9.0248e+00 -9.0248e+00  2e-07  2e-16  5e-08  4e-08
Optimal solution found.
person Dalek    schedule 01.08.2014
comment
Ницца! Вы можете подавить вывод cvxopt, поставив cvxopt.solvers.options['show_progress'] = False после импорта. Вот блокнот, в котором собрано все решение. - person Abraham D Flaxman; 02.08.2014