optimization

[Matt Hoffman]

Hi stan-dev,

I added a --point_estimate flag that runs Newton's method to find a
maximum of the joint log-likelihood, using finite-differences to
compute the Hessian and an eigendecomposition to adjust it to be
negative definite if necessary.

Is it all right if I push it? I know we don't really want to encourage
people to blindly fit point estimates, but sometimes it's useful to be
able to. For example, I'm currently interested in fitting a model by
maximum marginal likelihood where I can integrate out the hidden
variables.

Maybe we could leave it as an undocumented feature.

Matt

[Bob Carpenter]

Please push. (Being a Detroiter by birth, I'm down with the
notion of releasing and then doing a recall if the vehicle
turns out to be too dangerous on the road.)

Does it also use the same input, data, and output
file parameters? Presumably none of the other
existing params won't matter.

Are there any new params?

I'd like to add it to the doc -- it should be a whole
new part. We probably want to talk about integrating
out and maximum marginal likelihood, so if you have
a simple example winging about you can share, please
push that, too.

- Bob

[Ben Goodrich]

On Thursday, November 29, 2012 3:42:21 PM UTC-5, Bob Carpenter wrote:

Please push.  (Being a Detroiter by birth, I'm down with the
notion of releasing and then doing a recall if the vehicle
turns out to be too dangerous on the road.)

I agree it should be pushed, but I don't actually see what is uniquely dangerous about it. People maximize likelihoods all the time, and if not in Stan then they will do it some other way. Am I correct in assuming that it does the gradient by auto-diff and the Hessian by finite-diffing the auto-diffed gradient? If so, then that is a good contribution over R's optimization functions for example, which either take a function that produces an analytic gradient or do numeric gradients.

Ben

[Matt Hoffman]

It's not really that dangerous, except that it could lead to more
people asking basic questions on the mailing list about the
differences between point estimation and Bayesian inference. Then
again, that's arguably an opportunity to educate, which is arguably a
good thing.

You are correct that it's finite-diffing the auto-diff gradients. I'm
being pretty conservative and using a 4th-order finite difference
formula rather than a tiny perturbation size, so hopefully the
numerical problems won't get too hairy.

I'll wrap up the output format and push. Probably this afternoon.

Matt

> --
>
>

[Matt Hoffman]

Ok, I've pushed my additions. The algorithm should be reasonably
robust, but I haven't tested it incredibly thoroughly. Relevant flags
work pretty much as before, except that many become irrelevant.
--samples specifies the output file, which is a bit of a misnomer.
--save_warmup specifies whether to save the path from the initial
state to the mode. --point_estimate makes the thing go.

It never ceases to amaze me how well Newton's method works (when you
can apply it). When deciding what to pass to std::setw() when
outputting iteration number, it was really nice to be confident that
there would probably never be a circumstance where the number of
iterations would go above 99.

Matt

[Jiqiang Guo]

Do we want Rstan to support it?

Jiqiang

--

[Andrew Gelman]

I agree that if it's in the program it should be in the manual.

To Ben:  the point is that mle can be ok for a model with just a few parameter but not for big hierarchical models of the sort that people are likely to put in stan.  If we had VB capabilities, that would be cool…

--
 
 

[Bob Carpenter]

On 11/29/12 6:46 PM, Jiqiang Guo wrote:
> Do we want Rstan to support it?

That'd be nice, but it's totally up to you when
to do it.

Also, Daniel and I can probably help out with the
R integration so you don't have to do just this
interface stuff. Maybe we could start by sitting
down and seeing how everything's put together and
how the builds work. It's also good to make sure
at least two people understand every aspect of
a product so when people leave, things don't go to
hell.

- Bob

[Matt Hoffman]

> If we had VB capabilities, that would be cool…

It's on my list.

[Ben Goodrich]

On Thursday, November 29, 2012 5:44:55 PM UTC-5, Andrew Gelman wrote:

I agree that if it's in the program it should be in the manual.

To Ben:  the point is that mle can be ok for a model with just a few parameter but not for big hierarchical models of the sort that people are likely to put in stan.  If we had VB capabilities, that would be cool…

I still don't see that as a reason not to put various optimization-related stuff into Stan. If people are trying to maximize a likelihood when they have an incidental parameter problem or whatever, that is not Stan's fault. But there is presumably a set of well-defined problems where analytic gradients are not known and numeric gradients are too slow or too imprecise but they can be solved with auto-diff gradients. I'm sure it would be easy enough to add L-BFGS-B and other quasi-Newton or trust region options if the numeric Hessian wasn't working well for a problem.

Ben

[Andrew Gelman]

Tha's fine, I just would like a sentence saying that joint optimization generally doesn't make sense for hierarchical models.

--
 
 

[Ben Goodrich]

On Thursday, November 29, 2012 6:46:07 PM UTC-5, Jiqiang Guo wrote:

Do we want Rstan to support it?

I would think so, eventually. But I would favor accessing it by a different function than stan() to reduce confusion, and having it return a plain list so that it was comparable to optim() or nlm().

Ben

[Matt Hoffman]

Agreed on all counts. There's no reason not to have an LBFGS option,
except that it's a bit more of a pain to implement correctly. But
obviously it's a much better option when there are lots of parameters.
And once we have it people could even use Stan to test and fit crazy
neural networks or big source separation models or whatever else
they're fitting point estimates for now.

Basically, I'm getting sick of writing model-specific inference code,
regardless of whether I'm using ML, VB, or MCMC. So I'll probably be
motivated to add these things as I need them.

Matt

> --
>
>

[Jiqiang Guo]

The good thing about a different function is that we don't need to differentiate results for sampling or just the optimized results (previously I have to deal with test_grad differently as there is no sampling).  

But the other problem is that for any case we need to compile the C++ code.  If we want to reuse the compiled code, we would want to feed a fitted results to do the optimization and vice versa.  Then it seems to be convenient to just use one function stan(). 

--

Jiqiang 

--
 
 

[Bob Carpenter]

It'd probably be a good idea to use a different
command in the command-line version of Stan, too.

I don't even know if anyone's using the command-line
stuff other than me for testing! We need user
surveys (just kidding).

- Bob

> --
>
>

[Bob Carpenter]

On 11/29/12 9:59 PM, Jiqiang Guo wrote:
> The good thing about a different function is that we don't need to differentiate results for sampling or just the
> optimized results (previously I have to deal with test_grad differently as there is no sampling).

Good point.

I just looked up the test_grad option, and it's another
one of the ... options. I see that it just prints, which
I think is fine -- there's no way from the API to dump out
the full vector data structure. I could do that if you think
it'd be useful or easier from the R side.

I still can't believe how much work this integration has
been -- it's brought up all sorts of issues I never thought
about. And you (Jiqiang) have added all sorts of nice
R-friendly touches.

> But the other problem is that for any case we need to compile the C++ code. If we want to reuse the compiled code, we
> would want to feed a fitted results to do the optimization and vice versa. Then it seems to be convenient to just use
> one function stan().

I hadn't thought about that either, but then I doubt lots of
people are going to be trying both MCMC and optimization
on the same problem, so I'd just go with whichever is easier
to write and maintain on the RStan side. Overall, I don't
think compilation time's that big of a big deal for the run sizes
people are doing.

- Bob

[Nigel/essence]

Not sure of the background of this discussion on optimization, but I see below talk of MLE and functions without analytic derivatives.

I do MLE optimisation in high dimensions, and currently use a heuristic algorithm (mainly DE) followed by BOBYQA. Neither needs derivatives. Once I have worked out the analytical derivative, I will use L-BFGS-B (I have converted from Fortran to Java) and maybe even a NUTS sampler, but still probably keep the initial heuristic algorithm.

My point is, if you are looking to do Newton type methods without derivatives, BOBYQA from one of optimizations gradfathers, Powell, is worth looking at. Originally Fortran, but there is a Java version from Apache.

[Matt Hoffman]

The nice thing about working in Stan is that you don't have to write
down (or even be able to write down) the derivatives—it uses automatic
differentiation to compute the derivatives of the log-joint function.
So if your model includes functions that are continuously
differentiable but whose derivatives don't have a closed form then
Stan can still compute those derivatives. It looks like BOBYQA is more
for the scenario where you only have an oracle that evaluates the
objective function; in Stan we have direct access to the objective
function, so we can take a more direct route.

[Nigel/essence]

OK, I think I was misunderstanding auto diff and finite diff. I assume auto diff does not involve any finite diffs.

Differentiation of MLE's can get quite messy, hence I currently use derivative-free methods. BOBYQA works remarkably well, but is not available in R AFAIK.:

see for example

http://www.mucm.ac.uk/Pages/Downloads/Technical%20Reports/11-03%20YA%20PC%20Parameter%20Estimation%20for%20Gaussian%20process%20emulators.pdf

[Andrew Gelman]

We do have an open question regarding numerical differential equation solvers, how to compute gradients in a way that isn't computationally killing.

> --
>

[Bob Carpenter]

The auto-diff does analytic derivatives on the program.
There's a pretty good description on the Wikipedia:

http://en.wikipedia.org/wiki/Automatic_differentiation

We use reverse mode with a templated C++ implementation of
what they call taping.

Basically, anything where you can code the probability
function in C++, we can differentiate it (as can the other
auto-diff packages). You can use loops, recursion, function
calls, whatever. You do have to either (a) use one of
the functions in our library, or (b) use something that
templates the variables for which you want to take the gradient.

The penalty in speed is about a factor of four over just
the probability function. This is independent of
dimension (unlike finite diffs, which is slow in every implementation and
even slower if you do something less naive than a single
point). Space usage can be bad for iterative
algorithms because it stores the entire expression tree leading
to the log probability in a forward pass before applying the
chain rule in a backward pass.

- Bob

On 11/30/12 2:29 PM, Nigel/essence wrote:
> OK, I think I was misunderstanding auto diff and finite diff. I assume auto diff does not involve any finite diffs.
>
> Differentiation of MLE's can get quite messy, hence I currently use derivative-free methods. BOBYQA works remarkably well, but is not available in R AFAIK.:
>
> see for example
>
> http://www.mucm.ac.uk/Pages/Downloads/Technical%20Reports/11-03%20YA%20PC%20Parameter%20Estimation%20for%20Gaussian%20process%20emulators.pdf
>
>
>
>
> On Friday, November 30, 2012 6:13:53 PM UTC, Matt Hoffman wrote:
>> The nice thing about working in Stan is that you don't have to write
>>

>> down (or even be able to write down) the derivatives�it uses automatic

[Nigel/essence]

Yes, I did automatic differentiation way back in 1980's, with user defined functions.

I'll bow out, but if you ever need to optimize a function without analytical derivatives, check out BOBYQA.

[Marcus Brubaker]

A worthwhile point here is that the definitive implementation of L-BFGS-B is licensed with a 3 clause BSD license here: http://users.eecs.northwestern.edu/~nocedal/lbfgsb.html

It would be a little work to interface with the fortran, but probably worthwhile.  In particular, the finite-diff to get the Hessian will be painfully slow and solving the resulting system may be infeasable with larger models.

Cheers,

Marcus

--

[Matt Hoffman]

That is indeed the definitive implementation. It's really annoying
that it's in Fortran. Not just because it's a pain, but more
importantly because I don't want to make Stan users install a Fortran
compiler.

There's no question that Newton's method is completely impractical in
many situations, and that we're (I'm) going to want LBFGS. It's not
that hard to implement the unconstrained version, it's just a little
bit of a pain. Fortunately we don't have to worry about box
constraints, because we're already transforming everything to an
unconstrained space.

Matt

> --
>
>

[Nigel/essence]

My understanding is that LBFGS (without bounds) had a relaxed license, if any, but LBFGS-B 2.0 had much more stringent copyright license, copyright to the publishers of the original journal, so you basically had to stay with the Fortran.

R uses LBFGS-B 2.0, last time I looked.

However, there is a LBFGS-B 3.0 which had a revised algorithm, released 2011, and this has the much more relaxed 3-clause license, and no copyright to any journal.

I have converted LBFGS-B 3.0 to Java, but it was not a trivial task. Problems with GOTO's and return values of arguments. Took me a few weeks of effort. So I would expect a similar conversion to C to take a few weeks.

The code for LBFGS is considerably different (and shorter) than LBFGS-B.