I'd like to use a Hessian based method for finding the set of parameters $\theta$ that maximizes function $f$ (it's actually a maximum likelihood problem). $f$ is nonconvex, and so the straightforward Newton-Raphson update method ($\theta_{k+1}=\theta_k - H^{-1}(\theta_k)\nabla{f}(\theta_k)$ ) generally fails to converge. This is because the method looks for any value of $\theta$ where the gradient is null, and so the iteration will equally look for maxima, minima or saddle points. Is there any general method to ensure that Newton-Raphson converges to a (local) maximum ?
I have given a try at homecooking one such method, and would love to hear from you if there is any consistent work that has looked at that.
The problem is that, for maximization, we always want the update to be in the direction of the gradient (not necessarily parallel to it, but in the direction of it). With NR, this will only occur if the function is locally concave, i.e. if all the eigenvalues of the hessian $H$ are negative, i.e. $H$ is negative definite.
My tweak is just to change $H^{-1}(\theta_k)$ in the update formula for a negative definite matrix. Using the eigenvector decomposition of $H$ ($Hv_k=\lambda_k v_k$, or $H=V\Lambda V'$), we change positive eigenvectors for negative ones, and leave the negative ones intact. We build a diagonal matrix $\Xi$, with diagonal elements
$$ \xi_k = - \alpha / \lambda_k \textrm{ with } 0\lt \alpha \leq 1 \textrm{ if } \lambda_k>0 $$
$$ \xi_k = \beta / \lambda_k \textrm{ with } 0\lt \beta \leq 1 \textrm{ if } \lambda_k \lt 0 $$
And then we simply replace the update by the modified rule : $\theta_{k+1}=\theta_k - (V \Xi V') \nabla f(\theta_k)$.
I used $\beta=1$ (or close to) to converge rapidly along the directions which display 'local concavity'. For the other directions I used a smaller step ($\alpha = 0.2$), since we need to carefully get away from the local 'convexity' and find the inflection point along that direction. Once we pass all inflection points and $H$ is locally concave, then the update is equivalent to a classical NR update (since then $\Xi=\beta \Lambda^{-1}$).
It does work for my problem, the iteration converges to a local maximum. Any comment on this ? Am I simply reinventing the wheel ?