# Adaptive gradient descent step size when you can't do a line search

I have an objective function $E$ dependent on a value $\phi(x, t = 1.0)$, where $\phi(x, t)$ is the solution to a PDE. I am optimizing $E$ by gradient descent on the initial condition of the PDE: $\phi(x, t = 0.0)$. That is, I update $\phi(x, t = 0.0)$ and then have to integrate the PDE to compute my residual. That means, if I were to do a line search for the gradient descent step size (call it $\alpha$), for every potential value of $\alpha$ I would have to integrate the PDE all over again.

In my case that would be prohibitively expensive. Is there another option for adaptive gradient descent step size?

I'm not just looking for mathematically principled schemes here (although of course that's better if something exists), but would be happy with anything that is generally better than a static step size.

Thanks!

• I don't think I want to modify the way I integrate the PDE at the moment, as for me that would be a major code rewrite. Also, it's not so much that the PDE is a tricky one, as that I have to solve it on a very dense grid in spacetime as I require very high numerical accuracy. – NLi10Me Jul 14 '16 at 18:41
• On the other hand, the BB method (which I wasn't familiar with) seems pretty good; all I have to do it keep track of the previous iteration's state and gradient and I get a second order approximation... that seems very nice. However, the derivation assumes a convex quadratic and my problem almost certainly isn't. Though, I am also certainly finding (and happy with) local rather than global minima. Do you know how well BB has performed on very high dimensional problems? – NLi10Me Jul 14 '16 at 18:43
• I guess what I meant about local minima is that, in the neighborhood of a local minimum, isn't any function approximately quadratic? I think my initial state $\phi^{(0)}(x, t = 0.0)$ is sufficiently close to a minimum, as for many instances I get smooth convergence even with the static step size. So, even though it's very high dimensional, and in general if you consider the whole search space the problem is non-convex/non-quadratic, could BB still be a good choice w/o line search? – NLi10Me Jul 14 '16 at 18:57
• The other "ingredients" to $E$ are experimental image data. $\phi(x, t = 1.0)$ attempts to warp one image to "match" the other (measured by some matching functional like L2 norm integrated over voxels). For some image pairs, I get smooth convergence with (my current choice of) static step size. For other image pairs, I get a lot of oscillating. The system has to be fully automated, so I can't go back and hand edit the step size for troublesome image pairs. – NLi10Me Jul 14 '16 at 19:02
• Right, I do have to solve the adjoint system to get the gradient (which is a nastier system and takes longer). Ok, I think I will try BB with backtracking line search. Thank you very much for the advice; my advisors are often hard to get ahold of and many of them aren't interested in the implementation so much as just the model. I'm finding the numerical methods are the crucial component to demonstrating whether a model is good or not in the first place, so thanks again I really appreciate it. – NLi10Me Jul 14 '16 at 19:15

I'll begin with a general remark: first-order information (i.e., using only gradients, which encode slope) can only give you directional information: It can tell you that the function value decreases in the search direction, but not for how long. To decide how far to go along the search direction, you need extra information (gradient descent with constant step lengths can fail even for convex quadratic problems). For this, you basically have two choices:

1. Use second-order information (which encodes curvature), for example by using Newton's method instead of gradient descent (for which you can always use step length $1$ sufficiently close to the minimizer).
2. Trial and error (by which of course I mean using a proper line search such as Armijo).

If, as you write, you don't have access to second derivatives, and evaluating the obejctive function is very expensive, your only hope is to compromise: use enough approximate second-order information to get a good candidate step length such that a line search needs only $\mathcal{O}(1)$ evaluations (i.e., at most a (small) constant multiple of the effort you need to evaluate your gradient).

One possibility is to use Barzilai--Borwein step lengths (see, e.g., Fletcher: On the Barzilai-Borwein method. Optimization and control with applications, 235–256, Appl. Optim., 96, Springer, New York, 2005). The idea is to use a finite difference approximation of the curvature along the search direction to get an estimate of the step size. Specifically, choose $\alpha_0>0$ arbitrary, set $g^0:=\nabla f(x^0)$ and then for $k=0,...$:

1. Set $s^k = -\alpha_k^{-1} g^k$ and $x^{k+1}=x^k+s^k$
2. Evaluate $g^{k+1}=\nabla f(x^{k+1})$ and set $y^k = g^{k+1}-g^{k}$
3. Set $\alpha_{k+1} = \frac{(y^k)^Ty^k}{(y^k)^Ts^k}$

This choice can be shown to converge (in practice very quickly) for quadratic functions, but the convergence is not monotone (i.e., the function value $f(x^{k+1})$ can be larger than $f(x^k)$, but only once in a while; see the plot on page 10 in Fletcher's paper). For non-quadratic functions, you need to combine this with a line search, which needs to be modified to deal with the non-monotonicity. One possibility is choosing $\sigma_k \in (0,\alpha_k^{-1})$ (e.g., by backtracking) such that $$f(x^k - \sigma_k g^k) \leq \max_{\max(k-M,1)\leq j\leq k} f(x^j) - \gamma \sigma_k (g^k)^Tg^k,$$ where $\gamma\in(0,1)$ is the typical Armijo parameter and $M$ controls the degree of monotonicity (e.g., $M=10$). There's also a variant that uses gradient values instead of function values, but in your case the gradient is even more expensive to evaluate than the function, so that doesn't make sense here. (Note: You can of course try to blindly accept the BB step lengths and trust your luck, but if you need any sort of robustness -- as you wrote in your comments -- that would be a really bad idea.)

An alternative (and, in my opinion, much better) approach would be to use this finite difference approximation already in the computation of the search direction; this is called a quasi-Newton method. The idea is to incrementally build an approximation of the Hessian $\nabla^2 f(x^k)$ by using differences of gradients. For example, you could take $H_0=\mathrm{Id}$ (the identity matrix) and for $k=0,\dots$ solve $$H_{k}s^{k} = -g^{k},\label{cc1}\tag{1}$$ and set $$H_{k+1} = H_k + \frac{(y^k-H_ks^k)^T(s^k)^T}{(s^k)^Ts^k}$$ with $y^k$ as above and $x^{k+1} = x^k +s^k$. (This is called Broyden update and is rarely used in practice; a better but slightly more complicated update is the BFGS update, for which -- and more information -- I refer to Nocedal and Wright's book Numerical Optimization.) The downside is that a) this would require solving a linear system in each step (but only of the size of the unknown which in your case is an initial condition, hence the effort should be dominated by solving PDEs to get the gradient; also, there exist update rules for approximations of the inverse Hessian, which only require computing a single matrix-vector product) and b) you still need a line search to guarantee convergence...

Luckily, in this context there exists an alternative approach that makes use of every function evaluation. The idea is that for $H_k$ symmetric and positive definite (which is guaranteed for the BFGS update), solving \eqref{cc1} is equivalent to minimizing the quadratic model $$q_k(s) = \frac12 s^T H_k s + s^T g^k.$$ In a trust region method, you would do so with the additional constraint that $\|s\| \leq \Delta_k$, where $\Delta_k$ is an appropriately chosen trust region radius (which plays the role of the step length $\sigma_k$). The key idea is now to choose this radius adaptively, based on the computed step. Specifically, you look at the ratio $$\rho_k := \frac{f(x^k)-f(x^k+s^k)}{f(x^k)-q_k(s^k)}$$ of the actual and predicted reduction in function value. If $\rho_k$ is very small, your model was bad, and you discard $s^k$ and try again with $\Delta_{k+1}<\Delta_k$. If $\rho_k$ is close to $1$, your model is good, and you set $x^{k+1}=x^k+s^k$ and increase $\Delta_{k+1}>\Delta_k$. Otherwise you just set $x^{k+1}=x^k+s^k$ and leave $\Delta_k$ alone. To compute the actual minimizer $s^k$ of $\min_{\|s\|\leq \Delta_k} q_k(s)$, there exist several strategies to avoid having to solve the full constrained optimization problem; my favorite is Steihaug's truncated CG method. For more details, I again refer to Nocedal and Wright.

• I'm just now looking at this again, and realize I have a question. In step three for the BB method you have $\alpha_{k+1} = \frac{(y^k)^Ty^k}{(y^k)^Ts^k}$; where $y^{k} = g^{k+1} - g^k$ and $s^k = -\alpha_k^{-1}g^k$. The numerator and denominator in the expression for $\alpha_{k+1}$ look like inner products. In my case, $g^k \in V^*$, where $V^*$ is a vector space with a non-trivial Riemannian metric: K. That is, $\langle g^k, g^k \rangle _{V^*} = \langle g^k, Kg^k \rangle_{L_2}$. Does that affect the definition of $\alpha_{k+1}$? – NLi10Me Aug 11 '16 at 23:14
• Yes, if you have a non-trivial vector space structure, you should respect that in the algorithms. In particular, you should distinguish between inner products of two functions in the same space (e.g., $y^k$ and $y^k$) and duality products between a function in the space and one in the dual space (e.g., $s^k$ and $y^k$) -- for the latter, you need to include the Riesz mapping to turn it into an inner product first. (This can be interpreted as preconditioning.) – Christian Clason Aug 12 '16 at 1:01
• Dr. Clason, I will be submitting a paper to ISBI 2017 detailing some experiments I have done using the BB + line search method for a diffeomorphic image registration task. Would you like to be included as an author on the manuscript? I have not written it yet, but I have most of the experiments either complete or underway. Please let me know. – NLi10Me Oct 1 '16 at 16:29
• @NLi10Me Thank you for the kind offer, but I haven't done anything that would merit coauthorship -- everything I wrote is standard textbook material. If you feel strongly about it, you can thank me for "helpful remarks about (whatever it is that helped)", but not even that would be required. Knowing that what I wrote was helpful is enough! – Christian Clason Oct 1 '16 at 17:17
• Sorry, you're right, that is a typo -- fixed! (The Armijo condition is often written as $f(x+\sigma s) - f(x) \leq \gamma\nabla f(x)^T(\sigma s)$, where $s$ is the search direction -- not necessarily the negative gradient -- and $\sigma$ the step size, which should make clearer what's going on.) – Christian Clason Oct 17 '16 at 19:09