# ADVI start with initialization

**URL:** <https://discourse.pymc.io/t/advi-start-with-initialization/359>\
**Category:** Questions\
**Created:** [September 21, 2017, 1:05pm UTC](https://discourse.pymc.io/t/advi-start-with-initialization/359 "2017-09-21T13:05:28Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![messcode](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/messcode/32/246_2.png) [@messcode](https://discourse.pymc.io/u/messcode)\
**Post date:** [September 21, 2017, 1:05pm UTC](https://discourse.pymc.io/t/advi-start-with-initialization/359/1 "2017-09-21T13:05:28Z")

</div>

Hey community!

I built a probabilistic model and I want to investigate the convergence property of latent parameter every N steps. The procedure can be summarized as follows:

```python
import pymc3 as pm
means = None
for dummy_i in xrange(10):
    with my_model:
        (means, std, elbos) = pm.variational.advi(n=N, start=means)
    check_means(means) 

```

I intend to track the model performance with iterations. However, I found it is relying on the step length N. In my case, when N=5000, I found that the model performance( AUC) is fluctuating at around 0.6. N =10000, the model performance is approaching 1 and the final model performance is around 0.9. It seems that stopping advi and continuing with the means of parameters gets a different result.  
So my problem is how can I stop iteration and continue?  
Thanks in advance.

---

<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 21, 2017, 1:27pm UTC](https://discourse.pymc.io/t/advi-start-with-initialization/359/2 "2017-09-21T13:27:38Z")

</div>

> [@messcode](#):
>
> So my problem is how can I stop iteration and continue?

The advi API you are using is a bit dated, if you update to master and use the new API, you can do the following:

```python
with pm.Model() as model:
    # define your model
    ... 
    # set up the advi inference
    inference = pm.ADVI()

# train your model
approx = inference.fit(n=10000)

```

Now you have an approximation class `approx`, which you can check the history of elbo, and also the meanfield approximation mu and sd

```python
elbos1 = -approx.hist # shape = (10000,)
# get the fitted value
gbij = approx.bij
cov = approx.cov.eval()
means = gbij.rmap(approx.mean.eval())
sds = gbij.rmap(np.diag(cov)**.5)

```

You can now check the fit by ploting the elbo or using the `means` and `sds`, if you want to further train the model, just do:

```python
inference.fit(n=20000)

```

and again do

```python
elbos1 = -approx.hist # shape = (10000,) + (20000,) = (30000,)
# get the fitted value
gbij = approx.bij
cov = approx.cov.eval()
means = gbij.rmap(approx.mean.eval())
sds = gbij.rmap(np.diag(cov)**.5)

```

to get the updated means an sds

---

<div class="post-metadata">

**Author:** ![messcode](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/messcode/32/246_2.png) [@messcode](https://discourse.pymc.io/u/messcode)\
**Post date:** [September 22, 2017, 6:06am UTC](https://discourse.pymc.io/t/advi-start-with-initialization/359/3 "2017-09-22T06:06:45Z")

</div>

Thanks for your detailed instruction @junpenglao . Very helpful! However I still have some problem.  
I defined my model by following codes:

```python
import pymc3 as pm
# y is a list of input matrices. and a, b alpha are parameters.
y = input_mats
a = args.a
b = args.b
alpha = 2
N = y[0].shape[0]
n_source = len(y)
dim = args.k
M = [y[c].shape[1] for c in range(n_source)]
prior_mu = np.zeros([N, dim])
thetas = {}
with pm.Model() as model:
    sigma = pm.InverseGamma('sigma', alpha=alpha, shape=n_source)
    x = pm.Laplace('x', mu=prior_mu, b=b, shape=(N, dim))
    for c in xrange(n_source):
        thetas['theta'+str(c)] = pm.Dirichlet('theta'+str(c), 
                                        a=(a / dim) * np.ones((M[c], dim)), shape=(M[c], dim))
    for c in xrange(n_source):
        pm.Normal('y'+str(c), mu=theano.tensor.dot(x, thetas['theta'+str(c)].T), 
                  sd=pm.math.sqrt(sigma[c]) * np.ones((N, M[c])), observed=y[c])
    inference = pm.ADVI()

```

I can use ADVI by the old interface:

```python
with model:
    means, std, elbos = pm.variational.advi(n=10000)

```

This works fine for me. However, when I use the new API:

```python
approx = inference.fit(n=10000)

```

It takes a very long time to start iteration.

 ![Selection_129](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/1X/88aca1084da0d3b11819605cdb4e81a6b4466c9f.png)  
Here you can see that inference.fit() takes much longer time (40 mins) to start iteration.  
ps: My OS is ubuntu 16.04 and I use anaconda python 2.7 and cuda-8.0.

---

<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 22, 2017, 6:56am UTC](https://discourse.pymc.io/t/advi-start-with-initialization/359/4 "2017-09-22T06:56:48Z")

</div>

The two for loop is likely inefficient. You should try to rewrite it.

There might be ways to check which part of the initialization is slow using some profiling, @ferrine?

Also, regarding to your original goal (i.e., checking the convergence property of latent parameter every N steps), you can [write your own callback function](http://docs.pymc.io/notebooks/variational_api_quickstart.html#Checking-convergence).

Alternatively, you can copy the old variational API and change the code in [L181-L201](https://github.com/pymc-devs/pymc3/blob/master/pymc3/variational/advi.py#L181-L201) to check the convergence internally (the old API will be removed soon).

---

<div class="post-metadata">

**Author:** ![ferrine](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ferrine/32/10_2.png) [@ferrine](https://discourse.pymc.io/u/ferrine)\
**Post date:** [September 22, 2017, 8:04am UTC](https://discourse.pymc.io/t/advi-start-with-initialization/359/5 "2017-09-22T08:04:11Z")

</div>

I suppose theano does graph optimizations and that takes a lot of time

---

<div class="post-metadata">

**Author:** ![messcode](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/messcode/32/246_2.png) [@messcode](https://discourse.pymc.io/u/messcode)\
**Post date:** [September 22, 2017, 8:37am UTC](https://discourse.pymc.io/t/advi-start-with-initialization/359/6 "2017-09-22T08:37:02Z")

</div>

Thanks for your advice. The parameter `n_source` is often a small number, e.g. 2 or 3.  
I profile those two line using `%%prun` in ipython.  
The results are as follows:

```python
%%prun 
with model:
    means, std, elbos = pm.variational.advi(n=10000)

```

The output file is like:

 ![profile_advi](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/1X/7550d1112abc00457b6f9589d1165cb9babdea67.png)  
New API:

```python
approx = inference.fit(n=10000)

```

The output is:

 ![profile_fit](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/1X/e5ee11ea7c0f222c69cb6fa4a2cb308f6b1eb494.png)  
It seems pickle.py and {repr} take the majority of running time.

---

<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 22, 2017, 9:06am UTC](https://discourse.pymc.io/t/advi-start-with-initialization/359/7 "2017-09-22T09:06:34Z")

</div>

Try to do `theta.append()` and indexing it with `theta[c]` instead of `thetas['theta'+str(c)]`

---

<div class="post-metadata">

**Author:** ![messcode](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/messcode/32/246_2.png) [@messcode](https://discourse.pymc.io/u/messcode)\
**Post date:** [September 23, 2017, 8:15am UTC](https://discourse.pymc.io/t/advi-start-with-initialization/359/8 "2017-09-23T08:15:56Z")

</div>

I tried, but it still doesn’t work. I modified the old variational API as you suggested. Thanks 🙂