Solved – Using L-BFGS-B optimization

optimizationr

This is a general question about how the L-BFGS-B optimization algorithm works.

I have encountered some strange likelihoods in a model I was running (which uses optim from R, and the L-BFGS-B algorithm). The function optimizes over a parameter, which is constrained to 0-1 and maximizes the likelihood (Minimizes the negative log-likelihood I believe is what it technically does)

As I said earlier, I found some odd results and since the model does not take long to run, I ran all possible values of the parameter and plotted the Likelihood for each value, the results looks like so:

enter image description here

The red line indicates the value that optim suggested, and the blue line is the max value I found by doing all possible values.

My question is then why is this happening?

Does the L-BFGS-B employ some kind of Newton-Raphson search algorithm, ie. it only finds local maxima / minima?

Or does the systematic process it goes through just not fit in a situation where the likelihood probability is as seen in the graph?

Thanks,

Sam

Best Answer

I will assume that the x-axis in your plot refers to indices of a vector of possible parameter values rather than the parameter values themselves (in which case as @mpiktas says there may be an issue with not respecting the constraint in your exhaustive search).

If this assumption is correct then the problem is that the optimal solution is on the boundary of the parameter space - i.e. from your exhaustive search the best value of the parameter is $1$. I suspect this may be compounded by your choice of starting value (is it $0$?).

I find it helpful in diagnosing issues with optim to include a print statement in the score function so I can follow which parameter values are being used at each function call. This may give you an intuition as to why it is not exploring certain regions of the parameter space.

Related Solutions

GLM – Does Log Likelihood in GLM Guarantee Convergence to Global Maxima?

The definition of exponential family is:

$$ p(x|\theta) = h(x)\exp(\theta^T\phi(x) - A(\theta)), $$

where $A(\theta)$ is the log partition function. Now one can prove that the following three things hold for 1D case (and they generalize to higher dimensions--you can look into properties of exponential families or log partition):

$ \frac{dA}{d\theta} = \mathbb{E}[\phi(x)]$
$ \frac{d^2A}{d\theta^2} = \mathbb{E}[\phi^2(x)] -\mathbb{E}[\phi(x)]^2 = {\rm var}(\phi(x)) $
$ \frac{ \partial ^2A}{\partial\theta_i\partial\theta_j} = \mathbb{E}[\phi_i(x)\phi_j(x)] -\mathbb{E}[\phi_i(x)] \mathbb{E}[\phi_j(x)] = {\rm cov}(\phi(x)) \Rightarrow \Delta^2A(\theta) = {\rm cov}(\phi(x))$

The above result prove that $A(\theta)$ is convex(as ${\rm cov}(\phi(x))$ is positive semidefinite). Now we take a look at likelihood function for MLE:

\begin{align} p(\mathcal{D}|\theta) &= \bigg[\prod_{i=1}^{N}{h(x_i)}\bigg]\ \exp\!\big(\theta^T[\sum_{i=1}^{N}\phi(x_i)] - NA(\theta)\big) \\ \log\!\big(p(\mathcal{D}|\theta)\big) &= \theta^T\bigg[\sum_{i=1}^{N}\phi(x_i)\bigg] - NA(\theta) \\ &= \theta^T[\phi(\mathcal{D})] - NA(\theta) \end{align}

Now $\theta^T[\phi(\mathcal{D})]$ is linear in theta and $-A(\theta)$ is concave. Therefore, there is a unique global maximum.

There is a generalized version called curved exponential family which would also be similar. But most of the proofs are in canonical form.

Solved – Maximum likelihood estimation for state space models using BFGS

I will give you an answer from my experience developing the R package stsm that is introduced in this document. By default ithe package uses the function optimize, which is a combination of a golden section search and successive parabolic interpolation. Other line search algorithms, such as the zoom line search algorithm will do the job similarly well. The algorithm in function optimize is convenient because it does not require the gradient of the function to be minimized. Here is some pseudo code to compute the optimal step size:

fcn.step <- function(step, model, direction)
{
  trial.pars <- model$pars + step * direction
  -logLik(model, pars = trial.pars)
}
s <- line.search(f = fcn.step, interval = c(0, 1))

fcn.step is the function to be minimized, the only parameter is the step size; the parameters of the model are fixed to the values from the last iteration of the BFGS optimization algorithm. The negative of the log-likelihood function is minimized with respect to the step size; direction is the optimal direction vector computed at the current iteration of the BFGS algorithm, $\tilde{\pi(\psi)}$ in your notation.

It is also important the choice of the upper limit of the interval where the optimal value of the step size s is searched. By default, the package stsm calculates the maximum value of s that is compatible with positive values for the variance parameters updated in equation (3). If this value is lower than $1$ then it is taken as the upper bound passed to the line search algorithm, otherwise this bound is set equal to $1$. You may look at function step.maxsize in the package or the function calc.constraint in this implementation of the BFGS algorithm.

With this approach you can apply the BFGS algorithm without the risk of reaching a local optimum with negative variance parameters. Otherwise, I would recommend you either the L-BFGS-B algorithm, which allows setting a lower bound equal to zero in the parameter space where the algorithm searches for the optimum; or a reparameterization of the model, e.g. variance = theta^2 or variance = exp(theta) and maximize the likelihood with respect to theta.

As regards the EM algorithm, as you mention, it is slow to converge. Strangely enough, it converges slower as it approaches the local optimum. For some insight into this issue and a modified EM algorithm you may look at this document. I am the author of this document, you may send me an e-mail if you have any questions about it or if you want a copy (the link is not always available).

Optimization of the likelihood function of state space models can be hard in practice. Some enhancements to the general-purpose optimization algorithm are recommended. For example, one of the parameters of the model can be concentrated out of the likelihood function; maximum likelihood in the frequency domain enjoys some advantages from a computational point of view and may provide good values to be used as starting values in the optimization of the time domain likelihood function.

Best Answer

Related Solutions

GLM – Does Log Likelihood in GLM Guarantee Convergence to Global Maxima?

Solved – Maximum likelihood estimation for state space models using BFGS

Related Question