# Some questions on the simplex transform

**URL:** <https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220>\
**Category:** Modeling\
**Tags:** performance\
**Created:** [July 12, 2017, 7:03pm UTC](https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220 "2017-07-12T19:03:22Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![aaronjg](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aaronjg/32/3768_2.png) [@aaronjg](https://discourse.mc-stan.org/u/aaronjg)\
**Post date:** [July 12, 2017, 7:03pm UTC](https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220/1 "2017-07-12T19:03:22Z")

</div>

I was looking at implementing the simplex type in a model, rather than the non identified softmax transform. I had some questions about the simplex transform and how it scales to large simplexes.

- The source code refers to a “centered stick breaking process.” Is this the same as the regular stick-breaking process referred to in the manual?
- Why is the transform done in the real scale rather than the log scale? It seems like log-scale might be more stable, especially if you are interested in `log(p)` in the end anyway. But maybe I’m overlooking something?
- Given the way the simplex is constructed, it seems like the ordering of the variables might be important for numerical stability/efficiency. Is it?
- It looks like with a simplex of size `k` would have the derivative propegate through `k-1` nodes, but if it were constructed through recursive splitting it would only have to propegate through `2*log(k)` nodes, which might be more numerically stable for large simplexes.

Is my intuition correct here, or are these wrong or unimportant?

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 13, 2017, 6:42pm UTC](https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220/2 "2017-07-13T18:42:26Z")

</div>

Yes, it’s the stick-breaking process in the model, but I implemented it in a centered way so that (0, 0, …, 0) on the unconstrained scale leads to a symmetric simplex (i.e., representing a uniform distribution).

All of our transforms are down to the natural constrained scales. It would be possible to provide a log\_simplex data structure either as a built-in or as a manual addition.

Yes, the ordering can matter when you have K in the thousands. There was a thread about this on discourse but I can’t find it.

You can work with completely unconstrained values (log odds) directly through something like our `categorical_logit` distribution. Why would you want to go to just the log scale?

Splitting in a binary rather than linear fashion sounds like it could be promising if the Jacobians are more stable. They could still be arranged so that everything’s lower-triangular and the determinant remains tractable. Usually underflow’s not an issue, though, if that’s what you’re worrying about.

It would be easy enough to try. The code for the transforms is all very modular in the implementations (look for `simplex_constrain` and `simplex_unconstrain` functions).

---

<div class="post-metadata">

**Author:** ![aaronjg](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aaronjg/32/3768_2.png) [@aaronjg](https://discourse.mc-stan.org/u/aaronjg)\
**Post date:** [July 14, 2017, 5:55pm UTC](https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220/3 "2017-07-14T17:55:33Z")

</div>

I wasn’t aware of the `categorical_logit`, and I was looking for it’s equivalent for `multinomial.`

Though, the categorical logit seems to be operating on unconstrained `k` rather than `k-1` scale.  
`CategoricalLogit(y|β) = Categorical(y|softmax(β))` - which is valuable for the specification of my model now, but less so when I move to the `k-1` parameterization to resolve the identifiability issues.

I’ll keep the binary tree thing in mind if underflow becomes an issue. I presume this would start to show up as divergent transitions if it becomes an issue?

It actually might be interesting to allow arbitrary trees in the simplex construction, and then the unconstrained parameters could be meaningful. This sort of stuff comes up as “balances” in compositional data in geology and biology.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 17, 2017, 10:21pm UTC](https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220/4 "2017-07-17T22:21:38Z")

</div>

Yes, `categorical_logit` takes `K`-vectors, so it’s not identified. It’d be easy enough to add a `K - 1` version.

We do run into problems out in the 1000s for simplexes and this balanced tree construction should be more effective there at making things symmetric. The trick in implementing it would be getting all the boundary conditions right in the recursion.

---

<div class="post-metadata">

**Author:** ![aaronjg](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aaronjg/32/3768_2.png) [@aaronjg](https://discourse.mc-stan.org/u/aaronjg)\
**Post date:** [July 18, 2017, 4:40pm UTC](https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220/5 "2017-07-18T16:40:57Z")

</div>

> [@Bob\_Carpenter](#):
>
> We do run into problems out in the 1000s for simplexes and this balanced tree construction

What sort of problems do you run into? I’m using a simplex with ~500-2000 parameters and want to keep an eye out for it. I am getting more divergent transitions than when I was using the weekly identifiable model but I’m not sure if that’s due to the parameter transformation or some other related change… It would be good to know what specifically to keep an eye out for - and if it turns out to be a problem, I can look into coding up the balanced tree type transformation.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 19, 2017, 10:08pm UTC](https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220/6 "2017-07-19T22:08:07Z")

</div>

See: [https://github.com/stan-dev/stan/issues/2273](https://github.com/stan-dev/stan/issues/2273)

---

<div class="post-metadata">

**Author:** ![aaronjg](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aaronjg/32/3768_2.png) [@aaronjg](https://discourse.mc-stan.org/u/aaronjg)\
**Post date:** [July 19, 2017, 10:31pm UTC](https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220/7 "2017-07-19T22:31:36Z")

</div>

Thanks! With working test case as well. Looks like Jussi was running into that with optimizing, does this come up with the HMC as well?

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 19, 2017, 10:46pm UTC](https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220/8 "2017-07-19T22:46:39Z")

</div>

I haven’t checked.

---

<div class="post-metadata">

**Author:** ![aaronjg](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aaronjg/32/3768_2.png) [@aaronjg](https://discourse.mc-stan.org/u/aaronjg)\
**Post date:** [July 20, 2017, 3:06am UTC](https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220/9 "2017-07-20T03:06:57Z")

</div>

Thanks. I’m running it now - it’s a slow test case - and will post back here when it finishes.

```
Gradient evaluation took 0.66874 seconds
1000 transitions using 10 leapfrog steps per transition would take 6687.4 seconds.
```

---

<div class="post-metadata">

**Author:** ![aaronjg](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aaronjg/32/3768_2.png) [@aaronjg](https://discourse.mc-stan.org/u/aaronjg)\
**Post date:** [July 21, 2017, 8:20am UTC](https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220/10 "2017-07-21T08:20:45Z")

</div>

I didn’t see any major issue running that model with HMC. The treedepth was very high at first but settled down after about 100 samples. There were no major changes between the beginning and end of the simplex in terms of estimated masses, Rhat, or parameter bias.

 ![](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/2bb524351e3f6013b8f0963b1391b75afcdb8ef5.png)

 ![](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/29dcd2fd8c134e9b5c2cddf0111beb5818a3eb10.png)

 ![](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/2bb524351e3f6013b8f0963b1391b75afcdb8ef5.png) ![](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/3e595cbe1ddd168e824372e7b5ccd7136cb3a9d7.png)

 ![](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/a812ec2fe5e512ec403bc594cb10678819236e19.png)

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 21, 2017, 11:44pm UTC](https://discourse.mc-stan.org/t/some-questions-on-the-simplex-transform/1220/11 "2017-07-21T23:44:19Z")

</div>

Good news and glad that you’re testing. I wonder why optimization is less robust. It may have something to do with the way L-BFGS approximates the Hessian interacting poorly with the simplex transform or with adaptation not being fast enough (I’d think it’d be the opposite and the problem would be too-fast adaptation).