# Any way to make Stan competitive with Tensorflow for maximum likelihood?

**URL:** <https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230>\
**Category:** Algorithms\
**Created:** [March 26, 2019, 11:03am UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230 "2019-03-26T11:03:32Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![jeffpollock9](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jeffpollock9/32/3826_2.png) [@jeffpollock9](https://discourse.mc-stan.org/u/jeffpollock9)\
**Post date:** [March 26, 2019, 11:03am UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/1 "2019-03-26T11:03:32Z")

</div>

I have been playing with some toy examples in stan and tensorflow and was wondering if there was a way I could speed up the stan version I have posted here:

> **[Maximum likelihood estimation with tensorflow probability and pystan](https://jeffpollock9.github.io/maximum-likelihood-estimation-with-tensorflow-probability-and-pystan/)**
>
> I'm quite excited about tensorflow 2 and have been using the alpha quite a lot recently. I'm also really enjoying a lot of the functionallity in tensorflow probability which I show a little of here. I installed both of these via: pip install...

It is a simple logistic regression with 1 million observations and 250 predictors which on my laptop tensorflow runs in about 7 seconds and stan about 100.

Any comments would be really appreciated!

---

<div class="post-metadata">

**Author:** ![rok\_cesnovar](https://avatars.discourse-cdn.com/v4/letter/r/7bcc69/32.png) [@rok\_cesnovar](https://discourse.mc-stan.org/u/rok_cesnovar)\
**Post date:** [March 26, 2019, 4:34pm UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/2 "2019-03-26T16:34:46Z")

</div>

What is the difference if you dont include the compilation of the Stan model?

---

<div class="post-metadata">

**Author:** ![jeffpollock9](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jeffpollock9/32/3826_2.png) [@jeffpollock9](https://discourse.mc-stan.org/u/jeffpollock9)\
**Post date:** [March 26, 2019, 5:17pm UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/3 "2019-03-26T17:17:57Z")

</div>

Hi @rok_cesnovar. The timing already does not include the compilation of the stan model, the compilation is done in:

```python
stan_model = pystan.StanModel(stan_file, extra_compile_args=["-O3", "-march=native"])

```

I just double checked by running `htop` and looking for the `cc1plus` process.

I am only timing the `optimizing` part:

```python
start = tm.time()
stan_mle = stan_model.optimizing(stan_data, init="0")
end = tm.time()

```

---

<div class="post-metadata">

**Author:** ![Charles\_Driver](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charles_driver/32/6524_2.png) [@Charles\_Driver](https://discourse.mc-stan.org/u/Charles_Driver)\
**Post date:** [March 26, 2019, 5:20pm UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/4 "2019-03-26T17:20:10Z")

</div>

I found the stan bfgs optimizer good in certain situations, but running large models with many parameters recently I changed to ‘ucminf’ (r optimizer) and cut the convergence time to 1/10th or so. Using a hacked up not yet robust stochastic gradient descent approach has reduced it still further.

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 27, 2019, 8:10am UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/5 "2019-03-27T08:10:57Z")

</div>

> [@jeffpollock9](#):
>
> It is a simple logistic regression with 1 million observations and 250 predictors which on my laptop tensorflow runs in about 7 seconds and stan about 100.

Can you check how many function evaluations each optimization makes?

> [@Charles\_Driver](#):
>
> I found the stan bfgs optimizer good in certain situations, but running large models with many parameters recently I changed to ‘ucminf’ (r optimizer) and cut the convergence time to 1/10th or so.

Since `ucminf` also uses BFGS ([CRAN - Package ucminf](https://cran.r-project.org/web/packages/ucminf/index.html)), do you know where the difference comes from? Different line search or just different stopping rule? Stan BFGS is using not that good defaults for different tolerances, so you may get faster results by changing those.

---

<div class="post-metadata">

**Author:** ![Charles\_Driver](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charles_driver/32/6524_2.png) [@Charles\_Driver](https://discourse.mc-stan.org/u/Charles_Driver)\
**Post date:** [March 27, 2019, 8:58am UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/6 "2019-03-27T08:58:02Z")

</div>

I have experimented with control parameters somewhat, though hardly a robust comparison. It’s not just a tolerance thing though. Ucminf combines bfgs and trust region approaches, it’s not a ‘regular’ bfgs - it was much better (on the limited number of higher dim problems i tried) than other r bfgs optimizers.

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 27, 2019, 9:18am UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/7 "2019-03-27T09:18:37Z")

</div>

> [@Charles\_Driver](#):
>
> it’s not a ‘regular’ bfgs

I’m now picky about the wording: BFGS refers to the update rule of the inverse Hessian and thus it seems it’s regular BFGS with with a trust region type monitoring of the line search. It helps if we can recognize where the performance differences come from: is it the update rule (there are others, but I’m still assuming ucminf is using regular BFGS) or is it the line search (ucminf has a better line search and trust region approach is part of that, but there can be other modifications, too). My guess is that the implementation of ucminf works better especially in the initial part when the quadratic approximation is bad and ucminf uses less function evaluations in the line search at that time. In my experience Stan bfgs is using most of the function evaluations for the two first line searches. It would be interesting to see a plot of function values vs. iterations comparing these optimization algorithms. Stan is also using limited memory version, which may affect convergence speed for high-dimensional problems, and you can try also increasing the memory limit, but I guess that would not bee enough to explain 10 fold difference in the number of function evaluations.

---

<div class="post-metadata">

**Author:** ![Charles\_Driver](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charles_driver/32/6524_2.png) [@Charles\_Driver](https://discourse.mc-stan.org/u/Charles_Driver)\
**Post date:** [March 27, 2019, 9:46am UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/8 "2019-03-27T09:46:27Z")

</div>

Fair point re terminology, I’m a bit sloppy there :) and yeah I already had the stan bfgs memory limit out to 50 or 100, you’re right that it helped things. I was going to show the comparison you mentioned, but it seems that save\_iterations=TRUE does not do anything, in Rstan at least.

`stano <<- optimizing(sm, data=standata, init=0, algorithm = 'LBFGS', verbose=TRUE, save_iterations=TRUE, tol_rel_grad=0, tol_grad=0,tol_param=0, history_size=50,tol_rel_grad=0,sample_file='stano.txt')`

---

<div class="post-metadata">

**Author:** ![jeffpollock9](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jeffpollock9/32/3826_2.png) [@jeffpollock9](https://discourse.mc-stan.org/u/jeffpollock9)\
**Post date:** [March 27, 2019, 10:27am UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/9 "2019-03-27T10:27:33Z")

</div>

Here is some information about the number of iterations and function evaluations:

```python
>>> stan_mle = stan_model.optimizing(stan_data, init="0", verbose=True, refresh=1)
Initial log joint probability = -693147
    Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes 
       1 -153067 398.553 8135.91 0.001 0.001 2   
       2 -134589 2.75453 5320.96 0.3396 0.3396 3   
       3 -120848 5.50538 1094.32 1 1 4   
       4 -120278 0.922224 350.409 0.9502 0.9502 5   
       5 -120154 0.445196 310.868 1 1 6   
       6 -119800 1.65806 553.019 1 1 7   
       7 -118957 4.55923 1122.38 1 1 8   
       8 -116593 13.8279 2205.73 1 1 9   
       9 -110364 38.0724 4064.61 1 1 10   
      10 -96659.6 95.6122 8271.67 1 1 11   
      11 -72972.9 125.326 14710.6 1 1 12   
      12 -48860.9 78.3897 20054.9 0.3134 1 15   
      13 -34399.5 12.6584 8131.33 0.1609 0.2987 17   
      14 -29751.8 10.3422 4594.59 0.1606 0.1606 18   
      15 -29040.7 1.035 1773.42 0.3283 1 20   
      16 -28746.1 1.40649 1296.17 0.3311 0.3311 21   
      17 -28163.3 5.90758 472.261 1 1 22   
      18 -28136.4 0.614142 286.503 0.2512 1 24   
      19 -28125.8 0.373414 224.702 0.2534 0.2534 25   
      20 -28107.3 1.53199 129.288 1 1 26   
      21 -28106.4 0.0614243 11.4379 0.1723 1 28   
      22 -28106.4 0.000643128 9.62061 0.174 0.174 29   
      23 -28106.4 0.00281231 0.10144 1 1 30   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance

```

When I was doing this I noticed that the first iteration takes a lot longer, could this be due to copying all the data? If you discount that, stan would be a lot faster than 100 seconds, however I’m not sure how to avoid it.

With the tensorflow implementation:

```python
>>> mle.num_objective_evaluations
<tf.Tensor: id=25332, shape=(), dtype=int32, numpy=106>
>>> mle.num_iterations
<tf.Tensor: id=25331, shape=(), dtype=int32, numpy=27>

```

I can see also that (as you said) the arguments/defaults are not quite the same in the two implementations, you can see [stan](https://pystan.readthedocs.io/en/latest/_modules/pystan/model.html#StanModel.optimizing) and [tensorflow](https://github.com/tensorflow/probability/blob/9a317db2e346ed9542fde64558fac8c6a5bd5f5e/tensorflow_probability/python/optimizer/lbfgs.py#L80).

Is there a way to expose the stan objective function and autodiff gradients to R or python? That way I could trial other optimizers like ucminf very easily. A quick google search didn’t reveal anything, but I think it would be a really useful thing to have.

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 27, 2019, 11:21am UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/10 "2019-03-27T11:21:17Z")

</div>

> [@jeffpollock9](#):
>
> Here is some information about the number of iterations and function evaluations:

Thanks. It seems that tensorflow is computing 3 times more evaluations, but much faster which is not surprise as it’s optimized for this kind of tasks. The difference could be just due to multithreading. There are pull requests for some speedups for Stan, too (in addition of GPU and multithreading).

> [@jeffpollock9](#):
>
> When I was doing this I noticed that the first iteration takes a lot longer, could this be due to copying all the data?

But not more than 50s?

> [@jeffpollock9](#):
>
> Is there a way to expose the stan objective function and autodiff gradients to R or python?

Yes. You can find an example for R at [log\_prob and grad\_log\_prob functions — log\_prob-methods • rstan](https://mc-stan.org/rstan/reference/stanfit-method-logprob.html) and for Python at [https://github.com/MichaelRiis/Python-ADVI/blob/master/Eightschool.ipynb](https://github.com/MichaelRiis/Python-ADVI/blob/master/Eightschool.ipynb)

---

<div class="post-metadata">

**Author:** ![jeffpollock9](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jeffpollock9/32/3826_2.png) [@jeffpollock9](https://discourse.mc-stan.org/u/jeffpollock9)\
**Post date:** [March 27, 2019, 11:59am UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/11 "2019-03-27T11:59:28Z")

</div>

@avehtari thanks for the links, I didn’t know about those rstan functions. I’ll try using those with ucminf when I have time and report back.

---

<div class="post-metadata">

**Author:** ![roualdes](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/roualdes/32/1594_2.png) [@roualdes](https://discourse.mc-stan.org/u/roualdes)\
**Post date:** [March 27, 2019, 4:12pm UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/12 "2019-03-27T16:12:32Z")

</div>

> [@jeffpollock9](#):
>
> tensorflow runs in about 7 seconds and stan about 100

These original numbers used tensorflow’s gpu based lbfgs method, correct?

When I run this tensorflow model on my laptop, with no gpu, it fits in about 19 seconds (23 iterations).

It’s my understanding that since the data were generated in tensorflow, then there’s no extra copying of data. Stan on the other hand is copying data around for this exercise.

In an attempt to better measure the time it took Stan to fit the model, without copying data around, I used CmdStan and inserted checkpoints like this one

`std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();`

I also matched tensorflow’s objective function tolerance settings, though not all of them because the two implementations don’t have all the same convergence criteria. I also matched history\_size, but this had less of an effect.

Under this setting Stan measures pretty consistently about 40 seconds (23 iterations).

Stan uses `double`s and tensorflow uses `float`s. Could this make up a 20 second difference?

---

<div class="post-metadata">

**Author:** ![jeffpollock9](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jeffpollock9/32/3826_2.png) [@jeffpollock9](https://discourse.mc-stan.org/u/jeffpollock9)\
**Post date:** [March 27, 2019, 4:28pm UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/13 "2019-03-27T16:28:47Z")

</div>

Hi @roualdes.

Yes I am using tensorflow with a gpu (nvidia 1060, not particularly fast).

I did generate the data in tensorflow on the gpu, but it gets copied to the cpu as a numpy array when I make the dictionary `stan_data` and I am not timing this part:

```python
stan_data = {"N": N, "P": P, "x": x.numpy(), "y": y[:, 0].numpy()}

```

I suspect that stan is doing a further copy of this data at the start of the `optimizing` call, or at least it is doing something which takes a long time before the iterations start moving smoothly (based on watching the output with `verbose=True` and `refresh=1`). I’m not sure if that is avoidable in pystan (or rstan).

For such big arrays I would expect float could be a lot faster (twice?) than double (can fit twice as many floats in a simd register and the cpu cache) but I don’t think it is possible to do that with stan right now. Not sure if there are any plans to allow different floating point arithmetic in stan but in some cases it could be quite useful.

You can choose various different floating point types in tensorflow, float (`tf.float32`) is by far the fastest on my gpu and is the default in tensorflow, though. If I have time I’ll try with double and run on my CPU later.

---

<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:** [March 29, 2019, 6:22pm UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/14 "2019-03-29T18:22:11Z")

</div>

> [@jeffpollock9](#):
>
> When I was doing this I noticed that the first iteration takes a lot longer, could this be due to copying all the data?

Maybe in PyStan. Data transfer can be expensive compared to L-BFGS in problems this size if done poorly. We found that in just reading data into RStan when evaluating optimization problems at this scale. The I/O was dwarfing the time to fit in L-BFGS.

---

<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:** [March 29, 2019, 7:32pm UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/15 "2019-03-29T19:32:58Z")

</div>

I followed up on the Jeff Pollock’s [original blog post](https://jeffpollock9.github.io/maximum-likelihood-estimation-with-tensorflow-probability-and-pystan/). Running in RStan on my 2012 Macbook Pro, it takes about 60s. (I turned off the Hessian and sampling.)

I didn’t try running in CmdStan as that will almost certainly be I/O dominated just getting the data in (250 predictors for 1M items is roughly 1GB with double-precision arithmetic). At least that’s what I found last time I tried CmdStan.

---

<div class="post-metadata">

**Author:** ![jeffpollock9](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jeffpollock9/32/3826_2.png) [@jeffpollock9](https://discourse.mc-stan.org/u/jeffpollock9)\
**Post date:** [March 29, 2019, 10:59pm UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/16 "2019-03-29T22:59:30Z")

</div>

Hi @Bob_Carpenter. Thanks for all the comments - really helpful.

I’ve updated the post and ran the example using rstan - it runs in about 30 seconds on my laptop!

Looking forward to running this again with `bernoulli_logit_glm`, `map_rect`, and the gpu stuff once it is available!

---

<div class="post-metadata">

**Author:** ![stevebronder](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stevebronder/32/17242_2.png) [@stevebronder](https://discourse.mc-stan.org/u/stevebronder)\
**Post date:** [March 30, 2019, 7:36pm UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/17 "2019-03-30T19:36:25Z")

</div>

> [@jeffpollock9](#):
>
> You can choose various different floating point types in tensorflow, float ( `tf.float32` ) is by far the fastest on my gpu and is the default in tensorflow, though.

We talked about this in stan but it would be kind of a tedious rewrite

---

<div class="post-metadata">

**Author:** ![Erik\_Strumbelj](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/erik_strumbelj/32/1710_2.png) [@Erik\_Strumbelj](https://discourse.mc-stan.org/u/Erik_Strumbelj)\
**Post date:** [July 23, 2019, 6:26am UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/18 "2019-07-23T06:26:37Z")

</div>

> [@jeffpollock9](#):
>
> Looking forward to running this again with `bernoulli_logit_glm` , `map_rect` , and the gpu stuff once it is available!

You can try out the GPU-optimized bernoulli\_logit\_glm\_lpdf if you don’t mind using cmdstan with an experimental Stan Math:

> **[GitHub - bstatcomp/gpu-stan-paper-materials: Replication scripts, measurement...](https://github.com/bstatcomp/gpu-stan-paper-materials)**
>
> Replication scripts, measurement data, visualization scripts, and installation instructions for the paper "GPU-based Parallel Computation Support for Stan". - GitHub - bstatcomp/gpu-stan-...

> **[GPU-based Parallel Computation Support for Stan](https://arxiv.org/abs/1907.01063)**
>
> This paper details an extensible OpenCL framework that allows Stan to utilize heterogeneous compute devices. It includes GPU-optimized routines for the Cholesky decomposition, its derivative, other matrix algebra primitives and some commonly used...

Logistic regression is actually one of the toy examples in the paper. I’m very interested in how this compares to TensorFlow on a single GPU. Experiments show that Stan GLMs on the GPU are about 50x faster than not using the GLM primitives and running on the CPU.

@rok_cesnovar, @stevebronder, what’s the ETA on all this being released in Stan Math? The GPU stuff will at some point in the near future be available through rstan?

---

<div class="post-metadata">

**Author:** ![rok\_cesnovar](https://avatars.discourse-cdn.com/v4/letter/r/7bcc69/32.png) [@rok\_cesnovar](https://discourse.mc-stan.org/u/rok_cesnovar)\
**Post date:** [July 23, 2019, 8:01am UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/19 "2019-07-23T08:01:03Z")

</div>

I have been historically bad at predicting the ETA for all the GPU features so I will instead list what is left to be merged:

- add a double template to matrix\_cl (see [https://github.com/stan-dev/math/pull/1281](https://github.com/stan-dev/math/pull/1281), done, reviewed & approved, waiting due to Jenkins Windows issues)
- add the var template to matrix\_cl (I think Steve is close to finishing that branch, then I hope a week of iterating with reviews if we will be efficient)
- add caching to matrix\_cl with the double template type (the above two lines should make this step easier, the last caching branch was too confusing and difficult to review so Steve split this in three PRs)
- add GPU GLMs one by one (these are ready and waiting, but need caching merged)

2.21 should definitely have GPU GLMs, estimating when they get merged to develop is a bit harder to do. So 2.21 should have all the feature we promised at the last Stancon + GLMs and cov\_exp\_quad.

Rstan 2.19 already has GPU support for cholesky\_decompose (that was added for Stan Math 2.19).

---

<div class="post-metadata">

**Author:** ![stevebronder](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stevebronder/32/17242_2.png) [@stevebronder](https://discourse.mc-stan.org/u/stevebronder)\
**Post date:** [July 23, 2019, 9:35am UTC](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230/20 "2019-07-23T09:35:44Z")

</div>

> [@rok\_cesnovar](#):
>
> add the var template to matrix\_cl (I think Steve is close to finishing that branch, then I hope a week of iterating with reviews if we will be efficient)

This is pretty much done besides some of the copy functions

[Next page](https://discourse.mc-stan.org/t/any-way-to-make-stan-competitive-with-tensorflow-for-maximum-likelihood/8230.md?page=2)