Stochastic optimization problem, ValueError

Questions v5
1

Hi everyone,

I am trying to solve following problem:

image

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\model.py:2034, in Deterministic(name, var, model, dims, auto)
   1973 """Create a named deterministic variable.
   1974 
   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?

2

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

2 Likes