# Adding constraints with pm.Potential

**URL:** <https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639>\
**Category:** v5\
**Tags:** modeling\
**Created:** [January 16, 2024, 11:00am UTC](https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639 "2024-01-16T11:00:45Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![knoahcud](https://avatars.discourse-cdn.com/v4/letter/k/4af34b/32.png) [@knoahcud](https://discourse.pymc.io/u/knoahcud)\
**Post date:** [January 16, 2024, 11:00am UTC](https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639/1 "2024-01-16T11:00:45Z")

</div>

Hi,  
I’m quite new to pymc and want to add strict constraints to a linear regression model.

```python
import pymc as pm
import pandas as pd
import numpy as np
import pytensor.tensor as at

if __name__ == " __main__":

    df = pd.read_csv("BLR_lib_GU48-1234.csv", header=None, index_col=False)

    X = df.iloc[:, :-1]
    # Standardize the features
    stds = X.std().values
    X -= X.mean()
    X /= X.std()

    y_obs = df.iloc[:, -1]

    with pm.Model(coords={"predictors": X.columns.values}) as model:
        # Prior on error SD
        sigma = pm.HalfNormal("sigma", sigma=0.1)

        # coefficients
        betas = pm.Normal("beta", sigma=1, dims="predictors")
        expected_value = at.dot(X.values, betas) 

        # constraints
        con_1 = pm.math.ge(betas/stds[0] + betas/stds[1], 0.0)
        con_2 = pm.math.ge(betas[2]/stds[2], 0.0)
        con_3 = pm.math.le(betas[4]/stds[4], 0.0)
        con_4 = pm.math.ge(betas[2]/stds[2] + betas[3]/stds[3] + betas[4]/stds[4], 0.0)
        con_5 = pm.math.ge(betas[5]/stds[5], 0.0)
        con_6 = pm.math.le(betas[8]/stds[8], 0.0)
        con_7 = pm.math.ge(betas[5]/stds[5] + betas[6]/stds[6] + betas[7]/stds[7] + betas[8]/stds[8], 0.0)

        potential_1 = pm.Potential("con_1", pm.math.log(pm.math.switch(con_1, 1, 0.0)))
        potential_2 = pm.Potential("con_2", pm.math.log(pm.math.switch(con_2, 1, 0.0)))
        potential_3 = pm.Potential("con_3", pm.math.log(pm.math.switch(con_3, 1, 0.0)))
        potential_4 = pm.Potential("con_4", pm.math.log(pm.math.switch(con_4, 1, 0.0)))
        potential_5 = pm.Potential("con_5", pm.math.log(pm.math.switch(con_5, 1, 0.0)))
        potential_6 = pm.Potential("con_6", pm.math.log(pm.math.switch(con_6, 1, 0.0)))
        potential_7 = pm.Potential("con_7", pm.math.log(pm.math.switch(con_7, 1, 0.0)))

        y = pm.Normal('y', mu=expected_value, sigma=sigma, observed=y_obs.values)

        start = {"sigma": 1.0, "beta": np.zeros(len(X.columns))}
        idata = pm.sample(draws=2000, tune=1000, chains=4, step=[pm.NUTS(target_accept=0.99, max_treedepth=30)])

```

The regression is performed with standardized input data. The restrictions are specified for the unstandardized coefficients, which is why they need to be transformed. During execution, when strict constraints are applied (using `pm.math.log(pm.math.switch(constraint, 1, 0.0))`), I receive the error message:

```auto
pymc.exceptions.SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'sigma_log__': array(-2.67328751), 'beta': array([ 0.95591867, 0.33871937, 0.00713485, 0.73745527, 0.74412467,
       -0.82814385, 0.69393661, 0.62493708, 0.95804281, -0.80145671,
       -0.51079383, 0.11153434, -0.66721516])}

```

This error disappears when I use a low number for the false state in `pm.math.switch`. Does anyone know why this is happening and whether the constraints have been integrated into the model?

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [January 16, 2024, 11:18pm UTC](https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639/2 "2024-01-16T23:18:39Z")

</div>

Welcome!

I cannot run your code because it relies on the external csv file. But I assume that the starting point[s] is[are] running afoul of your constraints. If the sampler cannot find valid points to begin its initialization, it may well give up with the error you’re seeing.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 17, 2024, 8:21am UTC](https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639/3 "2024-01-17T08:21:00Z")

</div>

Hard constraints like that are unlikely to play well with NUTS, even if it finds a valid initial point. It’s better if you can softly penalize them instead.

Some of those also seem avoidable altogether by using positive priors

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [January 17, 2024, 10:02am UTC](https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639/4 "2024-01-17T10:02:06Z")

</div>

`ZeroSumNormal` prior might be useful, combined with something like `pm.Potential('positive_penalty', -(betas[4] + betas[8])` to punish draws with positive values in the positions of the beta vector you want to be negative.

I don’t think dividing by the standard deviation is relevant here, because it’s always strictly positive.

---

<div class="post-metadata">

**Author:** ![knoahcud](https://avatars.discourse-cdn.com/v4/letter/k/4af34b/32.png) [@knoahcud](https://discourse.pymc.io/u/knoahcud)\
**Post date:** [January 18, 2024, 8:31am UTC](https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639/5 "2024-01-18T08:31:20Z")

</div>

Thank you for all of your advices!

> [@cluhmann](#):
>
> If the sampler cannot find valid points to begin its initialization, it may well give up with the error you’re seeing.

I have tried using the `initval` argument in `pm.sample` and `pm.Normal`. Both doesn’t seem to work. When I check for the inital point with `model.initial_point()` before `pm.sample` it shows me the set initial values. Once the `pm.sample()` starts, other values are used. What could be the problem here?

> [@ricardoV94](#):
>
> Hard constraints like that are unlikely to play well with NUTS, even if it finds a valid initial point. It’s better if you can softly penalize them instead.
> 
> Some of those also seem avoidable altogether by using positive priors

I changed some priors resulting which leaves only 3 constraints. Is there a threshold where you would say that using hard constraints is still acceptable with NUTS or should i just soft threshold the remaining 3? When soft thresholding those constraints, which small number would be sufficient to emulate a hard constraint?

> [@jessegrabowski](#):
>
> `ZeroSumNormal` prior might be useful, combined with something like `pm.Potential('positive_penalty', -(betas[4] + betas[8])` to punish draws with positive values in the positions of the beta vector you want to be negative.
> 
> I don’t think dividing by the standard deviation is relevant here, because it’s always strictly positive.

I don’t understand how `ZeroSumNormal` can be used here. The constraints are all less or greater equal. Can you please explain the basic idea a bit further.  
The constraints are given in the unstandardized formulation to ensure numerical stability in the resulting equation.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 18, 2024, 8:55am UTC](https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639/6 "2024-01-18T08:55:47Z")

</div>

> [@knoahcud](#):
>
> Is there a threshold where you would say that using hard constraints is still acceptable with NUTS or should i just soft threshold the remaining 3? When soft thresholding those constraints, which small number would be sufficient to emulate a hard constraint?

NUTS will never behave well with hard constraints if the posterior is close enough to “violating them”. By soft-contraints I actually mean a continuous (differentiable) function that increasingly penalizes the posterior as it violates the constraints more and more. The shape and strength of this penalty will be something you have to investigate as it depends on the model/data.

The typical example from ML is to penalize very large coefficients with something like `sum(beta ** 2) * penalty_multiplier`.

In a bayesian setting, the ideal is to define constraints as part of the prior. But that is sometimes not possible or too cumbersome, in which case a Potential is a nice tool to just throw at the problem.

---

<div class="post-metadata">

**Author:** ![knoahcud](https://avatars.discourse-cdn.com/v4/letter/k/4af34b/32.png) [@knoahcud](https://discourse.pymc.io/u/knoahcud)\
**Post date:** [January 18, 2024, 10:10am UTC](https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639/7 "2024-01-18T10:10:57Z")

</div>

> [@ricardoV94](#):
>
> NUTS will never behave well with hard constraints if the posterior is close enough to “violating them”. By soft-contraints I actually mean a continuous (differentiable) function that increasingly penalizes the posterior as it violates the constraints more and more. The shape and strength of this penalty will be something you have to investigate as it depends on the model/data.
> 
> The typical example from ML is to penalize very large coefficients with something like `sum(beta ** 2) * penalty_multiplier`.
> 
> In a bayesian setting, the ideal is to define constraints as part of the prior. But that is sometimes not possible or too cumbersome, in which case a Potential is a nice tool to just throw at the problem.

So for example `con_1 = pm.math.ge(betas/stds[0] + betas/stds[1], 0.0)` could be penalized with `pm.Potential("con_1", pm.math.switch(con_1, 0, -((betas/stds[0] + betas/stds[1])**2*penalty_multiplier)` instead of a fixed penalty as done beforehand.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 18, 2024, 10:34am UTC](https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639/8 "2024-01-18T10:34:58Z")

</div>

Something like that, but probably the negative of the expression, as the Potential is added to the model logp. You want the logp to be lower when the constraint is violated more

---

<div class="post-metadata">

**Author:** ![knoahcud](https://avatars.discourse-cdn.com/v4/letter/k/4af34b/32.png) [@knoahcud](https://discourse.pymc.io/u/knoahcud)\
**Post date:** [January 18, 2024, 12:03pm UTC](https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639/9 "2024-01-18T12:03:25Z")

</div>

Isn’t the negative in the switch of `pm.Potential("con_1", pm.math.switch(con_1, 0, -((betas/stds[0] + betas/stds[1])**2*penalty_multiplier)` enough or am I missing something?

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [January 18, 2024, 12:45pm UTC](https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639/10 "2024-01-18T12:45:30Z")

</div>

My logic with ZeroSum was that you want some of the coefficients to be negative and some to be positive. If they all sum to zero, satisfying one set of constraints (say the negative ones) would automatically push the others parameter values towards satisfying the other constraints.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 18, 2024, 12:48pm UTC](https://discourse.pymc.io/t/adding-constraints-with-pm-potential/13639/11 "2024-01-18T12:48:51Z")

</div>

Sorry I didn’t see the minus before. That should be enough