Stochastic Programming with MCMC

Question

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.

Answer 1

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)$:

A notebook with all the code for this is here.