# Only one chain running, three are stuck

**URL:** <https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676>\
**Category:** v5\
**Tags:** modeling, pytensor\
**Created:** [August 7, 2023, 10:27am UTC](https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676 "2023-08-07T10:27:45Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![M\_K-C](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/m_k-c/32/6150_2.png) [@M\_K-C](https://discourse.pymc.io/u/M_K-C)\
**Post date:** [August 7, 2023, 10:27am UTC](https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676/1 "2023-08-07T10:27:45Z")

</div>

Hi all,

this problem is a bit of a strange one: When running my current model with `pm.sample()` only one of the four changes is exploring the posterior, while the other three are stuck at their initial values and don’t change.  
**About the model:**  
I am trying to infer parameters from a multidimensional dataset (4x2500) using a model I have written in PyTensor. A custom PyTensor function is part of the model. There are multiple combinations of parameters that could explain the data so the posterior may be rather flat at regions.

Everything runs fine, if I set `chains=1`.  
Are there any things that could lead to this effect of only one chain running? Maybe I need to specify the parallel chains in the custom PyTensor function?

Any help or ideas are highly appreciated…

---

<div class="post-metadata">

**Author:** ![M\_K-C](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/m_k-c/32/6150_2.png) [@M\_K-C](https://discourse.pymc.io/u/M_K-C)\
**Post date:** [August 14, 2023, 1:58pm UTC](https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676/2 "2023-08-14T13:58:41Z")

</div>

Update: This issue disappears when I use pm.Metropolis(). Therefore I think it is related to one of two things:  
a) the ‘flatness’ of the posterior that makes the NUTS sampler struggle somehow  
or  
b) an error in the grad() function inside my custom PyTensor function, as that one won’t be used by the Metropolis-sampler.

For now I’ll stick with Metropolis, but if you have any other suggestions, I’d be happy to hear them.

---

<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:** [August 14, 2023, 2:33pm UTC](https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676/3 "2023-08-14T14:33:54Z")

</div>

Is this issue unique to this model? Do other models sample well with multiple chains?

---

<div class="post-metadata">

**Author:** ![M\_K-C](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/m_k-c/32/6150_2.png) [@M\_K-C](https://discourse.pymc.io/u/M_K-C)\
**Post date:** [August 14, 2023, 2:43pm UTC](https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676/4 "2023-08-14T14:43:30Z")

</div>

Yes, it seems to be unique for this model with the custom pyTensor function.  
Other models run fine with multiple chains.

---

<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:** [August 14, 2023, 4:39pm UTC](https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676/5 "2023-08-14T16:39:54Z")

</div>

> [@M\_K-C](#):
>
> A custom PyTensor function is part of the model.

Can you expand on that? Maybe share some code?

---

<div class="post-metadata">

**Author:** ![M\_K-C](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/m_k-c/32/6150_2.png) [@M\_K-C](https://discourse.pymc.io/u/M_K-C)\
**Post date:** [August 15, 2023, 1:56pm UTC](https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676/6 "2023-08-15T13:56:16Z")

</div>

Yes, this function is the result of [another Discourse post](https://discourse.pymc.io/t/custom-pytensor-op-for-root-finding-with-multiple-inputs-and-outputs/11450/11):

```auto
class RootFinder(Op):
        
    def __init__ (self, t):
        
        self.t= t
        
    def make_node(self, D, a, b):

        outputs = [at.vector(dtype='float64')]
        return Apply(self, [D, a, b], outputs)
    
    def perform(self, node, inputs, outputs_storage):
        D, a, b = inputs
        
        outputs_storage[0][0] = root_finder_loop_perform(D, a, b, self.t)
    
    def grad(self, inputs, output_gradients):
        D, a, b = inputs
        x_list= self(D, a, b)
        
        x_grad_list = a*D/((a**2*b)/4 - D**2*x_list**2) + 2*a*D**3*x_list**2/((a**2*b)/4 - D **2*x_list** 2)**2
        D_grad_list = 4*a*(a **2*b*x_list+ 4*D** 2*x_list**3)/(a**2*b-4*D **2*x_list** 2)**2
        a_grad_list = -4*D*x_list*(a **2*b + 4*D** 2*x_list**2)/(a**2*b-4*D **2*x_list** 2)**2
        b_grad_list = -4*a**3*D*x_list/(a**2*b - 4*D **2*x_list** 2)**2
        
        grad_D = at.dot((-D_grad_list/x_grad_list), output_gradients[0])
        grad_a = at.dot((-a_grad_list/x_grad_list), output_gradients[0])
        grad_b = at.dot((-b_grad_list/x_grad_list), output_gradients[0])
       
        return grad_D, grad_a, grad_b

```

where, `root_finder_loop_perform()` is a function that returns the first n roots of a function. D, a and b are random scalar variables.The `grad()` function is expected to return three zero-dimensional outputs. I have confirmed that it works with grad\_verify() and it runs smoothly with a single chain.

I then call the function via:

```auto
rootfinder = RootFinder(t=t)
x_list = rootfinder(D, a, b)

```

and use the n-dimensional x\_list values in the next steps.  
As I said previously, this works fine for a single chain, but doesn’t seem to like multiple chains.  
Any help is highly appreciated!

---

<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:** [August 15, 2023, 3:03pm UTC](https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676/7 "2023-08-15T15:03:36Z")

</div>

What gets called in `root_finder_loop_perform`, and does it use multiple cores? There might be some multiprocessing scheduling issues resulting from each PyMC worker spawning lots of workers that don’t know about each other.

---

<div class="post-metadata">

**Author:** ![M\_K-C](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/m_k-c/32/6150_2.png) [@M\_K-C](https://discourse.pymc.io/u/M_K-C)\
**Post date:** [August 15, 2023, 3:44pm UTC](https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676/8 "2023-08-15T15:44:41Z")

</div>

It calls scipy.optimize.shgo() to find the n roots of a function. Is that known to use multiple cores?

---

<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:** [August 15, 2023, 5:20pm UTC](https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676/9 "2023-08-15T17:20:50Z")

</div>

[Very possibly](https://github.com/scipy/scipy/blob/c20c60f2a149c314611b02aff5652ad5665132b4/scipy/optimize/_shgo_lib/_vertex.py#L267)?

---

<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:** [August 15, 2023, 8:35pm UTC](https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676/10 "2023-08-15T20:35:33Z")

</div>

Does the problem persist when you set workers=1 on shgo? Might be barking up the wrong tree entirely, but that would be the test.

---

<div class="post-metadata">

**Author:** ![M\_K-C](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/m_k-c/32/6150_2.png) [@M\_K-C](https://discourse.pymc.io/u/M_K-C)\
**Post date:** [August 15, 2023, 9:00pm UTC](https://discourse.pymc.io/t/only-one-chain-running-three-are-stuck/12676/11 "2023-08-15T21:00:58Z")

</div>

Yes, I just tried that and it doesn’t solve the issue. But maybe the shgo algorithm is less effective with workers = 1 and hence lead to more problems. Maybe using pm.Metropolis() is the best way in this case…