There's a fun algorithm I learnt about recently, gradient descent with Sobolev preconditioning (or just preconditioning in general), that I felt like sharing (even though I still feel like there are aspects I don't totally understand).
So, in gradient descent, we're trying to find the minimum (or at least a local minimum) of some differentiable function f : ℝ^n -> ℝ. We start with some initial guess, x_0, and iterate, x_{n+1} = x_n - t_n ∇f (x_n), where the step size t_n is some real number (which may depend on x_n and f, or may just be constant, different implementations of gradient descent do this differently). We want x_{n+1} to be as close to the minimum as possible.
Consider the simple case where f is a quadratic function, f(x) = (x-b)ᵀ A (x-b)/2 for symmetric positive-definite A. A C2 function sufficiently near a local minimum can be approximated in this form. Wlog, by changing the coordinates, we can assume that A is diagonal and b is zero. Then ∇f(x_n) = A x_n, and x_{n+1} = x_n - t_n A x_n. Setting t_n to 1/(A_ii) will guarantee that the ith coordinate of x_{n+1} is zero, which is the value of the correct answer after the coordinate transform, but this only works for one axis at a time. If the entries of A are very different from each other, a step size t_n that works well for a small one will overshoot for a large one, and vice-versa. There is in general no value of the step size that will give fast convergence. (You could pick a different axis to minimize on different steps, but once you add back in the complexity of an arbitrary function f, that becomes much harder.)
So if A has entries of very different sizes (is ill conditioned) (for non-diagonal A, the condition is that its eigenvalues are of very different sizes), convergence is poor. Why not do another change of coordinates so that A isn't ill-conditioned any more? Scale x on the ith axis (in the diagonal case) by 1/√A_ii. Then f(x) becomes simply xᵀx/2, or x·x/2, and convergence is immediate. However, this change of coordinates takes more care than the rotation required to make A diagonal. ∇_x f(Mx) is not generally equal to (∇f)(Mx), rather it is Mᵀ(∇f)(Mx), so applying this re-scaling will materially change the algorithm. If we're effectively applying gradient descent to find the vextor Mx such that f(M⁻¹Mx) is minimal, then the iteration formula becomes
Mx_{n+1} = Mx_n - t_n M⁻¹ᵀ ∇f(M⁻¹Mx)
x_{n+1} = x_n - t_n (MᵀM)⁻¹ ∇f(x),
the question is then how to pick a good M (or (MᵀM)⁻¹ directly, as that's the only place it actually affects the algorithm) so that convergence is fast.
The matrix A which gives the best local quadratic approximation to an arbitrary function f is the Hessian (the matrix of second derivatives) of f, which is usually something we can explicitly calculate. In the diagonal quadratic case described previously, the matrix M which gave immediate convergence was √A. As A is symmetric, (MᵀM)⁻¹ simplifies to A⁻¹, and we can apply the same approach in general using the Hessian, giving the iteration formula x_{n+1} = x_n - (∇∇f(x_n))⁻¹ ∇f(x_n). For a quadratic function, this converges exactly in one step, and for an initial guess x_0 close enough to the true minimum that f is well approximated there by a quadratic, it gives much better convergence than plain gradient descent. It's equivalent to solving ∇f(x) = 0 by Newton's method, but this gives some indication of where it would fail: it isn't always finding a minimum, sometimes it's finding a maximum or saddle point instead. Also, if any of the (x_n)s are at or near inflection points (which will be all over the place in a high-dimensional space), the division may cause huge errors or fail entirely.
Instead of setting the preconditioning matrix (MᵀM)⁻¹ to exactly the inverse of the Hessian, which gives perfect conditioning near the minimum but sometimes doesn't work at all, and may be computationally expensive, it can be useful to set it instead to some value, designed to fit the particular application, that merely gets the eigenvalues of the preconditioned Hessian to be of generally the same scale rather than exactly the same.
Suppose f = g∘h, where f is ill-conditioned for gradient descent but g is not. It would be possible to minimize g(y) easily, but there may not be any x such that h(x) = argmin(g), and even if there is, x may be difficult to find. Instead of running gradient descent on g though, we can run preconditioned gradient descent on f, using the decomposition into g∘h to inform the preconditioning, generalizing the earlier case where f = (f ∘M⁻¹) ∘ M was decomposed into well-conditioned f ∘M⁻¹ and linear M. Rather than taking h to be linear, just take it to have a linear approximation, i.e. a derivative. Then the preconditioning matrix becomes ((∇h(x_n))ᵀ∇h(x_n))⁻¹. Even if the Jacobian matrix ∇h(x_n) isn't square, its product with its transpose will be. And of course if (∇h(x_n))ᵀ∇h(x_n) is possibly singular, or hard to compute, then an approximation can be used just like with the Hessian of the whole function f.
A special case of this is Sobolev preconditioning. Suppose the argument to f is (a discretization of) a function on some manifold Ω, and f depends on the gradient of this function. Then set h to the discrete gradient operator (which is actually linear), and the elements x'ᵀ (∇h)ᵀ∇h x of the matrix become ∫_Ω ∇_z x(z) · ∇_z x'(z) dz = ∫_∂Ω x(z) ∇_z x'(z) · n dz - ∫_Ω ∇²_z x(z) x'(z) dz, where the integral comes from taking the dot product of two functions. If the boundary term can be neglected (either because the manifold has no boundary, or because of the boundary conditions), this implies that the preconditioning matrix is just (-∇²)⁻¹, where the Laplacian is over the manifold Ω so after discretization it becomes just an n×n matrix, where n is the number of sample points. It may be a singular matrix (depending on boundary conditions), since the Laplacian of a constant function is 0, but in practice there are ways to work around this.
There's quite a nice intuition for what the Sobolev preconditioner is doing to make gradient descent faster. If the function x is already somewhat optimized, but at some point z_0 in its domain there is a gradient in f requiring x to change (i.e. d/dt f(x + t δ(z-z_0)) is non-zero), then after one simple gradient descent step, x(z_0) will change, which will change the gradients of x around z_0, then since f depends on the gradient of x, the next gradient descent step will update these adjacent points to compensate for the change, then the next step will change the points one step farther out from z_0, and so on. With a Sobolev preconditioner, instead of just updating x at z_0, the values at all the surrounding points are immediately dragged alond with it, since the effect of (-∇²)⁻¹ is to spread the input out like an electric field spreads around a set of charges. This is done in a way that minimizes the total change in the gradients of x.
I heard of this from a video about a particular application using an extension of Sobolev preconditioning.
stormborns
tangentwithtanner
titsay
gummygunk
janainaart
snowghoul
ellievsbear
annaxmalina
