Stochastic Gradient Descent

This tutorial illustrates how to use the stochastic_gradient_descent solver and different DirectionUpdateRules in order to introduce the average or momentum variant, see Stochastic Gradient Descent.

Computationally, we look at a very simple but large scale problem, the Riemannian Center of Mass or Fréchet mean: for given points $p_i ∈\mathcal M$, $i=1,…,N$ this optimization problem reads

$$\operatorname*{arg\,min}_{x∈\mathcal M} \frac{1}{2}\sum_{i=1}^{N} \operatorname{d}^2_{\mathcal M}(x,p_i),$$

which of course can be (and is) solved by a gradient descent, see the introductionary tutorial or Statistics in Manifolds.jl. If $N$ is very large, evaluating the complete gradient might be quite expensive. A remedy is to evaluate only one of the terms at a time and choose a random order for these.

We first initialize the packages

using BenchmarkTools, Colors, Manopt, Manifolds, PlutoUI, Random

and we define some colors from Paul Tol

begin
    black = RGBA{Float64}(colorant"#000000")
    TolVibrantOrange = RGBA{Float64}(colorant"#EE7733") # Start
    TolVibrantBlue = RGBA{Float64}(colorant"#0077BB") # a path
    TolVibrantTeal = RGBA{Float64}(colorant"#009988") # points
end;

We next generate a (little) large(r) data set

begin
    n = 5000
    σ = π / 12
    M = Sphere(2)
    x = 1 / sqrt(2) * [1.0, 0.0, 1.0]
    Random.seed!(42)
    data = [exp(M, x, random_tangent(M, x, Val(:Gaussian), σ)) for i in 1:n]
    localpath = join(splitpath(@__FILE__)[1:(end - 1)], "/") # files folder
    image_prefix = localpath * "/stochastic_gradient_descent"
    @info image_prefix
    render_asy = false # on CI or when you do not have asymptote, this should be false
end

false

render_asy && asymptote_export_S2_signals(
    image_prefix * "/center_and_large_data.asy";
    points=[[x], data],
    colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal]),
    dot_sizes=[2.5, 1.0],
    camera_position=(1.0, 0.5, 0.5),
)

false

render_asy && render_asymptote(image_prefix * "/center_and_large_data.asy"; render=2);

PlutoUI.LocalResource(image_prefix * "/center_and_large_data.png")

Note that due to the construction of the points as zero mean tangent vectors, the mean should be very close to our initial point x.

In order to use the stochastic gradient, we now need a function that returns the vector of gradients. There are two ways to define it in Manopt.jl: either as a single function that returns a vector, or as a vector of functions.

The first variant is of course easier to define, but the second is more efficient when only evaluating one of the gradients.

For the mean, the gradient is

$$ gradF(x) = \sum_{i=1}^N \operatorname{grad}f_i(x) \quad \text{where} \operatorname{grad}f_i(x) = -\log_x p_i$$

which we define in Manopt.jl in two different ways: either as one function returning all gradients as a vector (see gradF), or – maybe more fitting for a large scale problem, as a vector of small gradient functions (see gradf)

F(M, x) = 1 / (2 * n) * sum(map(p -> distance(M, x, p)^2, data))

F (generic function with 1 method)

gradF(M, x) = [grad_distance(M, p, x) for p in data]

gradF (generic function with 1 method)

gradf = [(M, x) -> grad_distance(M, p, x) for p in data];

The calls are only slightly different, but notice that accessing the 2nd gradient element requires evaluating all logs in the first function. So while you can use both gradF and gradf in the following call, the second one is (much) faster:

x_opt1 = stochastic_gradient_descent(M, gradF, x)

3-element Vector{Float64}:
 0.7071067811865475
 0.0
 0.7071067811865475

@benchmark stochastic_gradient_descent($M, $gradF, $x)

BenchmarkTools.Trial: 4647 samples with 1 evaluation.
 Range (min … max):  848.750 μs …   8.276 ms  ┊ GC (min … max): 0.00% … 86.38%
 Time  (median):     984.558 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):     1.071 ms ± 619.767 μs  ┊ GC (mean ± σ):  7.06% ± 10.09%

  ▃█▄▁                                                          ▁
  ████▁▃▃█▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▆▆ █
  849 μs        Histogram: log(frequency) by time        5.9 ms <

 Memory estimate: 861.45 KiB, allocs estimate: 10030.

x_opt2 = stochastic_gradient_descent(M, gradf, x)

3-element Vector{Float64}:
 0.7071067811865475
 0.0
 0.7071067811865475

@benchmark stochastic_gradient_descent($M, $gradf, $x)

BenchmarkTools.Trial: 10000 samples with 5 evaluations.
 Range (min … max):  5.080 μs … 553.814 μs  ┊ GC (min … max):  0.00% … 94.64%
 Time  (median):     6.540 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   8.621 μs ±  24.630 μs  ┊ GC (mean ± σ):  16.53% ±  5.75%

   ▃▇█▇▅▃                                                 ▁▁▁ ▂
  ▇███████▆▇▇█▇▇▆▆▃▁▁▃▄▁▁▁▁▁▁▁▁▁▃▁▁▁▁▃▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▅▇███ █
  5.08 μs      Histogram: log(frequency) by time      28.9 μs <

 Memory estimate: 40.98 KiB, allocs estimate: 23.

x_opt2

3-element Vector{Float64}:
 0.7071067811865475
 0.0
 0.7071067811865475

This result is reasonably close. But we can improve it by using a DirectionUpdateRule, namely:

On the one hand MomentumGradient, which requires both the manifold and the initial value, in order to keep track of the iterate and parallel transport the last direction to the current iterate. You can also set a vector_transport_method, if ParallelTransport() is not available on your manifold. Here, we simply do

x_opt3 = stochastic_gradient_descent(
    M, gradf, x; direction=MomentumGradient(M, x, StochasticGradient(zero_vector(M, x)))
)

3-element Vector{Float64}:
 0.7071067811865475
 0.0
 0.7071067811865475

MG = MomentumGradient(M, x, StochasticGradient(zero_vector(M, x)));

@benchmark stochastic_gradient_descent($M, $gradf, $x; direction=$MG)

BenchmarkTools.Trial: 10000 samples with 6 evaluations.
 Range (min … max):  4.950 μs … 467.745 μs  ┊ GC (min … max):  0.00% … 96.88%
 Time  (median):     6.484 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   8.247 μs ±  22.024 μs  ┊ GC (mean ± σ):  16.87% ±  6.23%

  ▁▄▆██▆▅▂                                                    ▂
  █████████████▇▆▅▄▄▁▁▃▃▁▁▄▄▁▁▁▁▁▁▃▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▆▆▇▇▇▅ █
  4.95 μs      Histogram: log(frequency) by time      27.3 μs <

 Memory estimate: 40.91 KiB, allocs estimate: 22.

And on the other hand the AverageGradient computes an average of the last n gradients, i.e.

x_opt4 = stochastic_gradient_descent(
    M, gradf, x; direction=AverageGradient(M, x, 10, StochasticGradient(zero_vector(M, x)))
);

AG = AverageGradient(M, x, 10, StochasticGradient(zero_vector(M, x)));

@benchmark stochastic_gradient_descent($M, $gradf, $x; direction=$AG)

BenchmarkTools.Trial: 10000 samples with 6 evaluations.
 Range (min … max):  4.467 μs … 395.773 μs  ┊ GC (min … max):  0.00% … 96.88%
 Time  (median):     6.467 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   8.428 μs ±  21.914 μs  ┊ GC (mean ± σ):  16.29% ±  6.18%

   ▂▃▄▇█▇▅▃▁                                                  ▂
  ▇███████████▇█▇▆▆▅▄▄▃▃▁▁▁▃▃▁▃▃▁▁▁▁▁▃▁▁▃▁▁▁▃▁▃▄▁▁▁▄▃▆▇████▇█ █
  4.47 μs      Histogram: log(frequency) by time      25.8 μs <

 Memory estimate: 40.91 KiB, allocs estimate: 22.

Note that the default StoppingCriterion is a fixed number of iterations which helps the comparison here.

For both update rules we have to internally specify that we are still in the stochastic setting (since both rules can also be used with the IdentityUpdateRule within gradient_descent.

For this not-that-large-scale example we can of course also use a gradient descent with ArmijoLinesearch, but it will be a little slower usually

fullGradF(M, x) = sum(grad_distance(M, p, x) for p in data)

fullGradF (generic function with 1 method)

x_opt5 = gradient_descent(M, F, fullGradF, x; stepsize=ArmijoLinesearch())

3-element Vector{Float64}:
  0.7071067811865475
 -3.9257729890729064e-16
  0.7071067811865475

AL = ArmijoLinesearch();

@benchmark gradient_descent($M, $F, $fullGradF, $x; stepsize=$AL)

BenchmarkTools.Trial: 6 samples with 1 evaluation.
 Range (min … max):  901.960 ms … 931.168 ms  ┊ GC (min … max): 7.32% … 6.66%
 Time  (median):     906.397 ms               ┊ GC (median):    7.33%
 Time  (mean ± σ):   910.165 ms ±  10.870 ms  ┊ GC (mean ± σ):  7.16% ± 0.32%

  █ █    █  █         █                                       █  
  █▁█▁▁▁▁█▁▁█▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  902 ms           Histogram: frequency by time          931 ms <

 Memory estimate: 711.73 MiB, allocs estimate: 9038266.