# Improving speed with large dimensional problem

**URL:** <https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345>\
**Category:** Questions\
**Created:** [September 19, 2017, 2:58am UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345 "2017-09-19T02:58:40Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![researcher](https://avatars.discourse-cdn.com/v4/letter/r/edb3f5/32.png) [@researcher](https://discourse.pymc.io/u/researcher)\
**Post date:** [September 19, 2017, 2:58am UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/1 "2017-09-19T02:58:40Z")

</div>

I previously asked about improving model sampling speed fitting about 14K observations to about 200 groups.

I am fitting a hierarchical model with Negative Binomial(alpha, mu) distributed likelihood (observations underneath level\_2). The priors for the parameters alpha and mu (level\_2) are Gamma(hyper\_alpha\_mu, hyper\_alpha\_sd) and Gamma(hyper\_mu\_mu, hyper\_mu\_sd) and the priors for the two levels above (level\_1, level\_0) are weakly informative HalfNormal distributions.

I am now fitting approximately 120K observations with 2000 classes at level\_2 and 7 classes at level\_1. The observations are count data on the order of magnitude of 10^5 - 10^7.

On a modestly powered laptop, the model goes through about 40 iterations a second during the initialization and it takes 3 hours to achieve 3000 samples using NUTS. When I try the sampling on a more powerful AWS instance, the sampling would take even longer.

Is there anything I can change in the model or the data to improve sampling speed?

See below for model specification.

Here is the dataset.

[anonymized\_dataset.csv](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/1X/0e5c50b8b18de7de68d1f625f7660e10e99af78a.csv) (3.0 MB)

```
import pandas as pd
import pymc3 as pm

def create_association(df, parent_col, child_col):
    renamed = df.rename(
        columns={parent_col: 'parent', child_col: 'child'}
    )
    association = renamed[['parent', 'child']].drop_duplicates()
    association.set_index('child', inplace=True)
    association.sort_index(inplace=True)

    return association

data = pd.read_csv('anonymized_dataset.csv')

associations = {
    'level_0': {'parent_col': 'level_0', 'child_col': 'level_1'},
    'level_1': {'parent_col': 'level_1', 'child_col': 'level_2'},
    'level_2': {'parent_col': 'level_2', 'child_col': 'observation'}
}

for level in associations:
    associations[level] = create_association(data, **associations[level])

print('number of groups at level_0: ', len(associations['level_0']['parent'].unique()))
print('number of associations of level_1 to level_0: ', len(associations['level_0']['parent'].values))
print('number of groups at level_1: ', len(associations['level_1']['parent'].unique()))
print('number of associations of level_2 to level_1: ', len(associations['level_1']['parent'].values))
print('number of groups at level_2: ', len(associations['level_2']['parent'].unique()))
print('number of observations associated to level_2: ', len(associations['level_2']['parent'].values))

observed = pd.merge(
    associations['level_2'], data[['observation', 'counts']],
    left_index=True, right_on='observation'
)['counts']

with pm.Model() as model:
    hyper_hyper_alpha_mu = pm.HalfNormal(
        'hyper_hyper_alpha_mu',
        sd=10,
        shape=len(associations['level_0']['parent'].unique())
    )
    hyper_hyper_alpha_sd = pm.HalfNormal(
        'hyper_hyper_alpha_sd',
        sd=10,
        shape=len(associations['level_0']['parent'].unique())
    )

    hyper_hyper_mu_mu = pm.HalfNormal(
        'hyper_hyper_mu_mu',
        sd=10,
        shape=len(associations['level_0']['parent'].unique())
    )
    hyper_hyper_mu_sd = pm.HalfNormal(
        'hyper_hyper_mu_sd',
        sd=10,
        shape=len(associations['level_0']['parent'].unique())
    )

    hyper_alpha_mu = pm.HalfNormal(
        'hyper_alpha_mu',
        hyper_hyper_alpha_mu[associations['level_0']['parent'].values],
        shape=len(associations['level_1']['parent'].unique())
    )
    hyper_alpha_sd = pm.HalfNormal(
        'hyper_alpha_sd',
        hyper_hyper_alpha_sd[associations['level_0']['parent'].values],
        shape=len(associations['level_1']['parent'].unique())
    )
    hyper_mu_mu = pm.HalfNormal(
        'hyper_mu_mu',
        hyper_hyper_mu_mu[associations['level_0']['parent'].values],
        shape=len(associations['level_1']['parent'].unique())
    )
    hyper_mu_sd = pm.HalfNormal(
        'hyper_mu_sd',
        hyper_hyper_mu_sd[associations['level_0']['parent'].values],
        shape=len(associations['level_1']['parent'].unique())
    )

    mu = pm.Gamma(
        'mu',
        mu=hyper_mu_mu[associations['level_1']['parent'].values],
        sd=hyper_mu_sd[associations['level_1']['parent'].values],
        shape=len(associations['level_2']['parent'].unique())
    )
    alpha = pm.Gamma(
        'alpha',
        mu=hyper_alpha_mu[associations['level_1']['parent'].values],
        sd=hyper_alpha_sd[associations['level_1']['parent'].values],
        shape=len(associations['level_2']['parent'].unique())
    )

    y_est = pm.NegativeBinomial(
        'y_est',
        mu=mu[associations['level_2']['parent'].values],
        alpha=alpha[associations['level_2']['parent'].values],
        observed=observed
    )

with model:
    hierarchical_trace = pm.sample(3000, progressbar=True)
```

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [September 19, 2017, 12:50pm UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/2 "2017-09-19T12:50:45Z")

</div>

The first thing I would try is changing the Gamma to a HalfNormal and use a non-centered parameterization.

Do you observed any speed up of NUTS before and after the tuning sample?

---

<div class="post-metadata">

**Author:** ![researcher](https://avatars.discourse-cdn.com/v4/letter/r/edb3f5/32.png) [@researcher](https://discourse.pymc.io/u/researcher)\
**Post date:** [September 19, 2017, 12:55pm UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/3 "2017-09-19T12:55:32Z")

</div>

By tuning sample, do you mean the initialization sample for ADVI? If so, for ADVI, it was about 30 it/s for the initialization and 1.5 s/it for the sampling from NUTS.

 ![05 AM](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/1X/ff54af0d515a8fc2a728d4a47aeb825d03261d4d.png)

Is there any literature on how to use a non-centered parameterization?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [September 19, 2017, 2:25pm UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/4 "2017-09-19T14:25:59Z")

</div>

By tuning I meant the speed of NUTS before and after 500 samples (ie, the default number of tuning).  
I am asking because if there is a speed difference, it indicates the initialization from ADVI is not good, PyMC3 ends up tuning the mass matrix in NUTS anyway, which sometimes could be quite long.

There are quite a few blog posts on parameterization in Hierarchical models:  
[http://twiecki.github.io/blog/2017/02/08/bayesian-hierchical-non-centered/](http://twiecki.github.io/blog/2017/02/08/bayesian-hierchical-non-centered/)

[http://docs.pymc.io/notebooks/Diagnosing\_biased\_Inference\_with\_Divergences.html](http://docs.pymc.io/notebooks/Diagnosing_biased_Inference_with_Divergences.html)

> **[Metamorphosis and the Multilevel Model](http://elevanth.org/blog/2017/09/07/metamorphosis-multilevel-model/)**
>
> \[latexpage\]Abstract: The same multilevel model can be written different ways. Algebra don’t mind. But the differences do matter to the computer. A terrible fact about applied statistical infe…

> **[Hamiltonian Monte Carlo for Hierarchical Models](https://arxiv.org/abs/1312.0906)**
>
> Hierarchical modeling provides a framework for modeling the complex
> interactions typical of problems in applied statistics. By capturing these
> relationships, however, hierarchical models also introduce distinctive
> pathologies that quickly limit the...

This discussion on Stan forum is pretty great too: [http://discourse.mc-stan.org/t/non-centered-parameterization-on-variance-parameter/584/12](http://discourse.mc-stan.org/t/non-centered-parameterization-on-variance-parameter/584/12)

---

<div class="post-metadata">

**Author:** ![researcher](https://avatars.discourse-cdn.com/v4/letter/r/edb3f5/32.png) [@researcher](https://discourse.pymc.io/u/researcher)\
**Post date:** [September 19, 2017, 2:41pm UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/5 "2017-09-19T14:41:22Z")

</div>

Hmm, I do recall that the first 500 samples appeared faster under the NUTS sampling. Should I specify more initial samples during ADVI?

Thank you for the blog posts on parameterization; I’ll check those out and see if they help.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [September 19, 2017, 3:17pm UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/6 "2017-09-19T15:17:32Z")

</div>

If you have the first 500 samples faster, it could be the case that the sampler fall into some kind of funnel. You should run multiple chain (e.g., `njobs=4`) and check the trace. Did you see any divergent warning?

---

<div class="post-metadata">

**Author:** ![researcher](https://avatars.discourse-cdn.com/v4/letter/r/edb3f5/32.png) [@researcher](https://discourse.pymc.io/u/researcher)\
**Post date:** [September 19, 2017, 3:19pm UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/7 "2017-09-19T15:19:27Z")

</div>

I did not see any divergent warning. I’ll retry with njobs=4, but doesn’t that default to the number of cores automatically which would be 4?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [September 19, 2017, 3:20pm UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/8 "2017-09-19T15:20:57Z")

</div>

Nope, you need to specify the number of chains. 🙂

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [September 19, 2017, 3:22pm UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/9 "2017-09-19T15:22:05Z")

</div>

FYI, have a look at this post for speed related diagnostic.

---

<div class="post-metadata">

**Author:** ![researcher](https://avatars.discourse-cdn.com/v4/letter/r/edb3f5/32.png) [@researcher](https://discourse.pymc.io/u/researcher)\
**Post date:** [September 19, 2017, 3:22pm UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/10 "2017-09-19T15:22:14Z")

</div>

Ok got it. Thank you for your advice. I’ll reply after trying the reparametrization and checking divergence.

---

<div class="post-metadata">

**Author:** ![researcher](https://avatars.discourse-cdn.com/v4/letter/r/edb3f5/32.png) [@researcher](https://discourse.pymc.io/u/researcher)\
**Post date:** [September 22, 2017, 4:26pm UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/11 "2017-09-22T16:26:01Z")

</div>

Hello, none of the links posted demonstrate how to reparametrize distributions outside of the location-scale family. Is it still possible to reparametrize the half-normal to achieve the same thing?

---

<div class="post-metadata">

**Author:** ![mcavallaro](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mcavallaro/32/187_2.png) [@mcavallaro](https://discourse.pymc.io/u/mcavallaro)\
**Post date:** [September 27, 2017, 4:23pm UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/12 "2017-09-27T16:23:48Z")

</div>

> [@junpenglao](#):
>
> FYI, have a look at this post for speed related diagnostic.

Which post? it seems that a link is missing 🙂

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [September 27, 2017, 5:32pm UTC](https://discourse.pymc.io/t/improving-speed-with-large-dimensional-problem/345/13 "2017-09-27T17:32:38Z")

</div>

Oops, here is the link [Frequently Asked Questions](https://discourse.pymc.io/t/frequently-asked-questions/74/4?u=junpenglao)