# Stochastic Programming with MCMC

While there are far superior methods for solving deterministic LP problems (e.g. interior point algorithms), can MCMC be used to solve their stochastic variants?

By stochastic, I mean, for example, with distributions on coefficients and or constraints, and getting a posterior probability on the solution representing:

• Uncertainty in the solution given the uncertainty in the coefficients/constraints
• Uncertainty in the solution considering that MCMC gives you an estimate

For example:

$$\begin{equation*} \begin{array}{rrl} \mathbf{x}^* = \underset{\mathbf{x}}{\text{arg}\;\text{min}} & \mathbb{E}\,(c_1x_1 -3x_2)\\ \mbox{s.t.} & -x_1 +x_2 & \le b_1 \\ & x_1, x_2 & \geq 0 \end{array} \end{equation*}$$

and where:

• $$c_1 \sim N(2,0.5)$$
• $$b_1 \sim N(0,3)$$

I am inclined to believe that we could use MRF potentials to represent the constraints and objective function, but not sure how to formulate a sampling problem that would help approximate a solution to the Eq. above.

• Can you provide a specific example of a problem you would like to solve? – Abraham D Flaxman Jul 3 '14 at 17:22

PyMC2 can be combined with the LP solver of your choice to solve stochastic LP problems like this one. Here is code to do it for this very simple case. I've left a note on how I would change this for a more complex LP.

c1 = pm.Normal('c1', mu=2, tau=.5**-2)
c2 = -3
b1 = pm.Normal('b1', mu=0, tau=3.**-2)

@pm.deterministic
def x(c1=c1, c2=c2, b1=b1):
# use an LP solver here for a complex problem
arg_min = np.empty(2)
min_val = np.inf
for x1,x2 in [[0,0], [0, b1], [-b1, 0]]: # there are only three possible extreme points,
if -x1 + x2 <= b1 and x1 >= 0 and x2 >= 0: # so check obj value at each valid one
val = c1*x1 + c2*x2
if val < min_val:
min_val = val
arg_min = [x1,x2]

return np.array(arg_min, dtype=float)


Look at the weird joint distribution for $(x_1, x_2)$:

• Thanks! This is very interesting. I was thinking of solving the problem instead by defining three potentials: one for the cost function and one for the constraints. The potential for the constraints would output p=0 (constraints not satisfied) or p=1 (constraints satisfied), while the potential for the cost function would return c1*x1 + c2*x2 + ... + cn*Xn. My understanding is that such approach would not require any LP solver. Your solution is very different, and it makes me wonder if I am wrong about the way I was thinking of modeling the problem. – Amelio Vazquez-Reina Jul 31 '14 at 2:25
• By the way, when I asked the question on GitHub, one of the contributors told me that using potentials would also be one way to proceed. I have to say that I like your solution a lot though, since I can imagine one could potentially wrap any solver with your approach! (i.e. it could be used to generalize any type of optimization problem, not just LP). – Amelio Vazquez-Reina Jul 31 '14 at 2:28
• Thanks Abraham for your answer and for uploading the details to a notebook. After looking into this more carefully, I am reaching the following conclusion: Your solution is useful for looking at the statistics of the optimal solution under variations of the model parameters, but it does not optimize the expected cost/utility given the uncertainty in the parameters of the formulation. This is an important clarification. I wonder if MCMC (e.g. PyMC) can be used to solve the latter on a relatively general case. – Amelio Vazquez-Reina Aug 5 '14 at 15:01
• Note that minimizing the expected cost (alternatively maximizing the expected utility) is a common optimization criterion in stochastic programming – Amelio Vazquez-Reina Aug 5 '14 at 15:05