Stochastic optimization problem, ValueError

Hi everyone,

I am trying to solve following problem:


Here is my code:

import numpy as np
import pymc as pm
from scipy.optimize import linprog

with pm.Model() as model:
    c1 = pm.Normal('c1', mu=2, tau=.5**-2)
    b1 = pm.Normal('b1', mu=0, tau=3)

    def lp_solve(c1=c1, b1=b1):
        obj = [c1.eval(), -3]

        lhs_ineq = [[-1, 1]]
        rhs_ineq = [b1.eval()]

        bnd = [(0, float("inf")),
               (0, float("inf"))]

        opt = linprog(c=obj, A_ub=lhs_ineq, b_ub=rhs_ineq, bounds=bnd, method="revised simplex")
        sol = np.array([opt['x'][0], opt['x'][1]], dtype=np.float32).flatten()
        return sol

    sols = pm.Deterministic(name="solutions", var=lp_solve(c1, b1))

    trace = pm.sample(20000)

I got following error:

ValueError                                Traceback (most recent call last)
Cell In [84], line 18
     15     sol = np.array([opt['x'][0], opt['x'][1]], dtype=np.float32).flatten()
     16     return sol
---> 18 sols = pm.Deterministic(name="solutions", var=lp_solve(c1, b1))
     20 trace = pm.sample(20000)

File ~\AppData\Local\Programs\Python\Python310\lib\site-packages\pymc\, in Deterministic(name, var, model, dims, auto)
   1973 """Create a named deterministic variable.
   1975 Deterministic nodes are only deterministic given all of their inputs, i.e.
   2031 var: var, with name attribute
   2032 """
   2033 model = modelcontext(model)
-> 2034 var = var.copy(model.name_for(name))
   2035 if auto:
   2036     model.auto_deterministics.append(var)

ValueError: order must be one of 'C', 'F', 'A', or 'K' (got 'solutions')

Any ideas?

This is happening because you are mixing numpy and aesara operations. pm.Deterministic expects an aesara tensor, but you are giving it a numpy array (np.ndarray.copy takes a string that defines the storage format of the array, while aesara.tensor.copy takes a string that defines the name of the symbol tensor object, hence the error).

If you want to use a scipy optimizer inside your model, you will need to wrap it into an aesara Op. Full docs on Ops are here, and an example inside a PyMC model can be found here.

Basically you need to tell it what datatypes and shapes to expect as inputs and outputs, then put your lp_solve code inside the perform method.

Note that doing this will not allow you to use NUTS, because your Op won’t have a gradient (unless you give it one).