Using Runge-Kutta integration to increase the speed and stability of gradient descent? For a gradient descent problem with $\mathbf{x}\in \mathbb{R}^N$ I can evaluate the gradient $\mathbf{\nabla}_\mathbf{x} \in \mathbb{R}^N$ that reduces the least squares error, $y$. However, simply updating the position using $\mathbf{x'} = \mathbf{x} + \mathbf{\nabla}_\mathbf{x}$ converges very slowly to the global minimum of the least squares error (which is also the global minimum of the gradient magnitude, where the gradient is zero). I tried simply scaling up the step, i.e.  $\mathbf{x'} = \mathbf{x} + h\mathbf{\nabla}_\mathbf{x}$, however while this dramatically improves convergence times in some cases, it can become unstable in others (particularly when some of the components of $\mathbf{\nabla}_\mathbf{x}$ are much larger than others -- scaling up all components of the gradient can cause the gradient descent method to "climb up the side of a canyon" rather than descending the canyon, and the system can either oscillate or explode).
I would like to use the 3rd order Runge-Kutta method to follow the curvature of the gradient space, so that I can take larger steps without the system blowing up. I have applied this to simulating mass-spring systems before (using Runge-Kutta integration to integrate acceleration to find velocity, and velocity to find position) -- however I can't figure out how to apply it to this gradient descent problem.
I think I have some fundamental misunderstanding about how the Runge-Kutta methods work. They requires a function $f=(x, y)$ to be defined, which I believe computes the gradient of the curve at $x$. However I don't understand why $y$ needs to be supplied to the function -- isn't $y$ a function of $x$?
Can Runge-Kutta even be applied to the gradient descent problem? It seems like there should be a way to adapt Runge-Kutta to gradient descent, since each update step $\mathbf{x'} = \mathbf{x} + \mathbf{\nabla}_\mathbf{x}$ is basically an integration step. Is the step size $h$ simply the magnitude of the gradient, i.e. $h_i = |{\mathbf{\nabla}_{\mathbf{x}_i}}|$ and $\mathbf{y}_i = {\mathbf{\nabla}_{\mathbf{x}_i}} / h_i$?
If Runge-Kutta is not applicable here, can somebody please suggest a robust and fast gradient descent algorithm to try?
Some more detail: in the case of this problem, the gradient surface is fairly smooth, and quite strongly convex (there are few if any local minima that are not global minima), but the error surface is less convex. In other words, sometimes gradient descent will continue walking down the gradient slope in the direction of the global minimum of gradient, and the least squares error will increase temporarily before decreasing to the global minimum of least squares error. (The gradient is not computed from the least squares error measure itself, but using a different method that directly identifies the locally-best least squares solution, which moves the system closer to the globally-optimal least squares solution.) The gradient is therefore more reliable for gradient descent than the slope of the least squares error surface.
 A: First, gradient descent and Runge-Kutta methods solve different problems.

*

*Gradient descent is a method to find an extremum of $f(\mathbf x)$ by solving $\mathbf g(\mathbf x) = \nabla f(\mathbf x) = 0$. Gradient descent simply does
$$
\mathbf x_{n+1} = \mathbf x_{n} + \alpha_n \mathbf g(\mathbf x_n)
$$
with $\alpha_n$ being fixed or chosen smartly.

*Runge-Kutta methods are used to solve ODEs, that is solve an initial value problem
$$
\mathbf x'(t) = \mathbf F(t, \mathbf x(t))\\
\mathbf x(0) = \mathbf x_0.
$$
The simplest RK method is Euler's method which has quite similar to GD form
$$
\mathbf x_{n+1} = \mathbf x_n + (t_{n+1} - t_n) \mathbf F(t_n, \mathbf x_n)
$$
In other words, GD may be treated as Euler's method applied to an ODE
$$
\mathbf x'(t) = \pm \mathbf g(\mathbf x)\\ \tag{*}
\mathbf x(0) = \mathbf x_0.
$$
I used $\pm$ since $\alpha_n$ may be positive or negative (depending whether you're searching for a minimum or a maximum). ODE's are usually solved forward in time, so $t_{n+1} - t_n$ is positive.
The solution you're searching is the steady state $\mathbf x(\infty)$ at which the left hand side (and, consequently the right side) becomes zero. The correct sign also ensures that $\mathbf x(t)$ really tends to the steady state and not away from it.
Further I will assume that the correct sign is $+$.
You may use higher order RK methods for (*) problem. For example, the midpoint rule
$$
\mathbf x_{n+1/2} = \mathbf x_{n} + \frac{\Delta t_n}{2} \mathbf g(\mathbf x_n)\\
\mathbf x_{n+1} = \mathbf x_{n} + \Delta t_n \mathbf g(\mathbf x_{n+1/2})
$$
Higher order RK methods are known to be more accurate than the Euler's method. That is the numerical trajectory (formed by $\mathbf x_n$ sequence) is much closer to the true trajectory $\mathbf x(t)$, which is the true solution of (*). Unfortunately, you don't need this property. In fact you don't care how close your $\mathbf x_n$ are to the true trajectory, instead you're interested in how close are your $\mathbf x_n$ to $\mathbf x(\infty)$.
It is attractive to choose $\Delta t_n$ large, so one faster approaches to the $t = \infty$. Unfortunately, it does not work that way, because all explicit methods for ODEs (and any RK method is one of them) have a stability condition that restricts the largest step $\Delta t$. In fact even choosing $\Delta t$ close to that bound won't work either since the method will be oscillating forward and backward (exactly like GD does). Choosing $\Delta t$ which maximizes the convergence speed is quite nontrivial.
Another disappointing fact is the stiffness phenomenon. You probably know that there are pathological functions $f(\mathbf x)$ for which GD converges very slowly. Usually it happens when Hessian matrix of $f$ is badly conditioned. For these cases the corresponding systems (*) are (infamously) known in numerical integration as stiff problems. For these problems all explicit methods perform roughly the same - the limit for $\Delta t$ and the convergence speed is believed to be practically the same.
Stiff problems are often solved by implicit methods. Those methods cannot be converted to a GD-like method, since they require solving a nonlinear problem for each iteration. And this problem is roughly equivalent to the minimization problem itself. For example implicit Euler method has the form
$$
\mathbf x_{n+1} = \mathbf x_{n} + \Delta t_n \mathbf g(\mathbf x_{n+1}).
$$
Separating known $\mathbf x_n$ and unknown $\mathbf x_{n+1}$ gives a nonlinear problem for $\mathbf x_{n+1}$
$$
\mathbf G(\mathbf x_{n+1}) \equiv \mathbf x_{n+1} - \Delta t_n \mathbf g(\mathbf x_{n+1}) = \mathbf x_{n}
$$
This problem is only slightly simpler to solve than the original $\mathbf g(\mathbf x) = 0$.
Summarizing all above: using more precise methods for (*) won't get you to the solution faster. Instead you may want to use conjugated gradients method or other methods that are specialized for minimization problems, possibly involving more information about the function.
