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.

  • $\begingroup$ Can you provide a specific example of a problem you would like to solve? $\endgroup$ Jul 3, 2014 at 17:22

1 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)

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

joint distribution

A notebook with all the code for this is here.

  • $\begingroup$ 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. $\endgroup$ Jul 31, 2014 at 2:25
  • $\begingroup$ 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). $\endgroup$ Jul 31, 2014 at 2:28
  • 1
    $\begingroup$ 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. $\endgroup$ Aug 5, 2014 at 15:01
  • 1
    $\begingroup$ Note that minimizing the expected cost (alternatively maximizing the expected utility) is a common optimization criterion in stochastic programming $\endgroup$ Aug 5, 2014 at 15:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.