Is it well known that some optimization problems are equivalent to time-stepping?

Question

Given a desired state $y_0$ and a regularization parameter $\beta \in \mathbb R$, consider the problem of finding a state $y$ and a control $u$ to minimize a functional \begin{equation} \frac{1}{2} \| y - y_0 \|^2 + \frac{\beta}{2} \| u \|^2 \end{equation} subject to the constraint \begin{equation} Ay = u. \end{equation} where for simplicity we can think of $ y, y_0, u \in \mathbb R^n $ and $ A \in \mathbb R^{n \times n} $.

Forming the Lagrangian, looking for stationary points, and eliminating the control $u$ we get the first-order conditions \begin{align*} A^T \lambda &= y_0 - y \\ Ay &= \frac{1}{\beta} \lambda \end{align*} Premultiplying by $A$ in the first equation and $A^T$ in the second, we can write the normal equations \begin{align} (I + \beta A A^T) \lambda &= \beta A y_0 \\ (I + \beta A^T A) y &= y_0 \end{align} We can interpret these as single steps of backward Euler approximations to the differential equations \begin{align} \frac{\partial \lambda}{\partial b} &= -A A^T \lambda + A y_0, \quad \lambda(0) = 0 \\ \frac{\partial y}{\partial b} &= -A^T A y, \quad y(0) = y_0 \end{align} with pseudotimestep $\beta$.

My question: Is this connection well known? Is it discussed in standard treatments of either timestepping or optimization? (To me, it seems to provide some kind of intuitive connection between them.)

The idea seems simple enough that it has to be well known, but neither searching the literature or talking to people has given me a good source where this is discussed. The closest I've found is a paper by O. Scherzer and J. Weichert (J. Math Imaging Vision 12 (2000) pp. 43-63) which states the connection in the first sentence of the abstract (!) but does not provide any references or explore the connection in any depth.

Ideally I am looking for a reference that not only states the connection but also explores some consequences (for example, one could imagine preconditioning an optimization problem with a cheap forward Euler step).

Answer 1

As Jed Brown mentioned, the connection between gradient descent in nonlinear optimization and time stepping of dynamical systems is rediscovered with some frequency (understandably, since it's a very satisfying connection to the mathematical mind since it links two seemingly different fields). However, it rarely turns out to be a useful connection, especially in the context you describe.

In inverse problems, people are interested in solving the (ill-posed) operator equation $F(u)=y^\delta$ with $y^\delta$ not in the range of $F$. (Your optimal control problem can be seen as one instance of it with $F=A^{-1}$ and $y^\delta = y_0$.) Several regularization strategies (such as Tikhonov or Landweber) can be interpreted as a single pseudo-time step of a certain class. The idea is then to use the interpretation of the regularization parameter as a step length to obtain some (adaptive, a posteriori) choice rules for the parameter -- a fundamental problem in inverse problems -- and possibly to make multiple pseudo-time steps to approach the true, unregularized solution (similarly to numerical continuation). This is sometimes called continuous regularization, and is usually discussed in the context of level set methods; see, for example, Chapter 6.1 of Kaltenbacher, Scherzer, Neubauer: Iterative Regularization Methods for Nonlinear Ill-Posed Problems (de Gruyter, 2008).

A second context this idea repeatedly crops up in is nonlinear optimization: If you look at a gradient descent step for $\min_x f(x)$, $$ x^{k+1} = x^k - \gamma_k \nabla f(x^k),$$ then you can interpret this as a forward Euler step for the dynamical system $$ \dot x(t) = -\nabla f(x(t)),\qquad x(0) = x^0.$$ As Jed Brown pointed out, this at first glance yields only the not very surprising observation that this method converges, provided the pseudo-time steps $\gamma_k$ are small enough. The interesting part comes when you look at the dynamical system and ask yourself what properties the continuous solution $x(t)$ of the so-called gradient flow has (or should have), independent of the gradient descent, and whether that might not lead to more appropriate time stepping (and hence optimization) methods than standard Euler. Some examples off the top of my head:

1. Is there a natural function space in which the gradient flow lives? If so, your gradient step should be taken from the same space (i.e., the discretization should be conforming). This leads, e.g., to computing Riesz representations of the gradient with respect to different inner products (sometimes called Sobolev gradients) and, in practice, to preconditioned iterations that converge much faster.

2. Maybe $x$ should belong not to a vector space, but to a manifold (e.g., symmetric positive definite matrices), or the gradient flow should conserve a certain norm of $x$. In this case, you could try to apply structure-preserving time-stepping schemes (e.g., involving a pull-back with respect to an appropriate Lie group or a geometric integrator).

3. If $f$ is not differentiable but convex, the forward Euler step corresponds to a subgradient descent method which can be very slow due to step size restrictions. On the other hand, an implicit Euler step corresponds to a proximal point method, for which no such restrictions apply (and which thus have become very popular in, e.g., image processing).

4. In a similar vein, such methods can be significantly accelerated by extrapolation steps. One way of motivating these is by observing that standard first-order methods suffer from having to make many small steps close to minimizers, because the gradient directions "oscillate" (think of the standard illustration for why conjugate gradients outperform steepest descent). To remedy this, one can "dampen" the iteration by not solving a first-order dynamical system, but a damped second-order system: $$a_1 \ddot x(t) + a_2 \dot x(t) = -\nabla f(x(t))$$ for suitably chosen $a_1,a_2$. With the proper discretization, this leads to an iteration (known as Polyak's heavy ball method) of the form $$x^{k+1} = x^k - \gamma_k \nabla f(x^k) + \alpha_k (x^k - x^{k-1})$$ (with $\gamma_k,\alpha_k$ depending on $a_1,a_2$). Similar ideas exist for proximal point methods, see, e.g., the paper http://arxiv.org/pdf/1403.3522.pdf by Dirk Lorenz and Thomas Pock.

---

(I should add that to my knowledge, in most of these cases the interpretation as a dynamical system was not strictly necessary for the derivation or the convergence proof of the algorithm; one could argue that ideas like "implicit vs. explicit" or Lie derivatives are actually more fundamental than either dynamical systems or gradient descent methods. Still, it never hurts to have another view point to look at a problem from.)

---

EDIT: I just stumbled across an excellent example from the second context, where the ODE interpretation is used to deduce properties of Nesterov's extragradient method and suggest improvements: http://arxiv.org/pdf/1503.01243.pdf (Note that this is also an example of Jed Brown's point, in that the authors essentially rediscover the point 4 above without apparently being aware of Polyak's algorithm.)

EDIT 2: And as an indication how far you can take this, see page 5 of http://arxiv.org/pdf/1509.03616v1.pdf.

Answer 2

While I haven't seen the exact formulation that you have written down here, I keep seeing talks in which people "rediscover" a connection to integrating some transient system, and proceed to write down an algorithm that is algebraically-equilavent to one form or another of an existing gradient descent or Newton-like method, and fail to cite anyone else. I think it's not very useful because the conclusion is basically that "as long as you take small enough steps, the method eventually converges to a local minimum". Well, 2014 marks the 45th anniversary of Philip Wolfe's paper showing how to do this in a principled way. There is also good theory for obtaining q-quadratic or q-superlinear convergence from pseudotransient continuation and related methods like Levenberg-Marquardt.

If you want an instance of this rediscovery using a Newton-like formulation for solving algebraic equations (i.e., classical pseudotransient continuation) from a mathematician with more than 600 papers (so maybe he'll prove things you find interesting), look at the "Dynamical Systems Method" by A.G. Ramm [1].

If the intuition gained by considering a transient system led to practical algorithms that were either faster or more reliable, I think we'd see highly-cited articles on that subject. I think it's no mystery that Nocedal and Wright has over 13000 citations while Ramm's book has about 80 (mostly self-citations).

[1] I can advise you not to inform Prof. Ramm that his DSM is algebraically-equivalent to something that has been in countless engineering packages for decades or you may get yourself yelled out of the room. #gradstudentmemories

Answer 3

If ODE methods can contribute to optimization, is there a really simple example problem to show this ?
A straw man: is there an ODE solver that does a reasonable job on
$\qquad \dot{ x } = - \nabla f( x ) $
or $\quad \ddot{ x } = \beta \dot{ x } - \alpha \nabla f( x ) \ \ $ as Christian Clason suggests
for $f$ say the Rosenbrock function, in 2d or 10d ? If that's silly, does anyone have a better straw man ?
(Note "reasonable", not "competitive with state-of-the-art optimizers". I imagine one needs decreasing step sizes / tolerance, and maybe a stiff solver.)

In practice, "too big" steps are much more problematic than "too small" — oscillations are messy.
I'd have thought naively that control theory could help. Numerical Recipes p. 915 describes
PI adaptive stepsize control for ODEs, but I don't know if this is used in practice.