Illustration of the Gradient of a Second Order Difference

This example explains how to compute the gradient of the second order difference mid point model using adjoint_Jacobi_fields.

This example also illustrates the PowerManifold manifold as well as ArmijoLinesearch.

We first initialize the manifold

exportFolder = joinpath(
    @__DIR__, "..", "..", "docs", "src", "assets", "images", "tutorials"
using Manopt, Manifolds

and we define some colors from Paul Tol

using Colors
black = RGBA{Float64}(colorant"#000000")
TolVibrantBlue = RGBA{Float64}(colorant"#0077BB") # points
TolVibrantOrange = RGBA{Float64}(colorant"#EE7733") # results
TolVibrantCyan = RGBA{Float64}(colorant"#33BBEE") # vectors
TolVibrantTeal = RGBA{Float64}(colorant"#009988") # geo

Assume we have two points $x,y$ on the equator of the Sphere $\mathcal M = \mathbb S^2$ and a point $y$ near the north pole

M = Sphere(2)
p = [1.0, 0.0, 0.0]
q = [0.0, 1.0, 0.0]
c = mid_point(M, p, q)
r = shortest_geodesic(M, [0.0, 0.0, 1.0], c, 0.1)
[c, r]

Now the second order absolute difference can be stated as (see [Bačák, Bergmann, Steidl, Weinmann, 2016])

\[d_2(x,y,z) := \min_{c ∈ \mathcal C_{x,z}} d_{\mathcal M}(c,y),\qquad x,y,z∈\mathcal M,\]

where $\mathcal C_{x,z}$ is the set of all mid points $g(\frac{1}{2};x,z)$, where $g$ is a (not necessarily minimizing) geodesic connecting $x$ and $z$.

For illustration we further define the point opposite of

c2 = -c

and draw the geodesic connecting $y$ and the nearest mid point $c$, namely

T = [0:0.1:1.0...]
geoPts_yc = shortest_geodesic(M, r, c, T)
nothing #hide

looks as follows using the asymptote_export_S2_signals export

    render = asyResolution,
    curves = [ geoPts_yc ],
    points = [ [x,y,z], [c,c2] ],
    colors=Dict(:curves => [TolVibrantTeal], :points => [black, TolVibrantBlue]),
    dotSize = 3.5, lineWidth = 0.75, cameraPosition = (1.2,1.,.5)
render_asymptote("SecondOrderData.asy"; render=2)

Three points \$p,r,q\$ and the midpoint \$c=c(p,q)\$ (blue)

Since we moved $r$ 10% along the geodesic from the north pole to $c$, the distance to $c$ is $\frac{9\pi}{20}\approx 1.4137$, and this is also what

costTV2(M, (p, r, q))

returns, see costTV2 for reference. But also its gradient can be easily computed since it is just a distance with respect to $y$ and a concatenation of a geodesic, where the start or end point is the argument, respectively, with a distance. Hence the adjoint differentials adjoint_differential_geodesic_startpoint and adjoint_differential_geodesic_endpoint can be employed, see ∇TV2 for details. we obtain

(Xp, Xr, Xq) = ∇TV2(M, (p, r, q))

When we aim to minimize this, we look at the negative gradient, i.e. we can draw this as

   points = [ [x,y,z], [c,c2] ],
   colors=Dict(:tvectors => [TolVibrantCyan], :points => [black, TolVibrantBlue]),
   dotSize = 3.5, lineWidth = 0.75, cameraPosition = (1.2,1.,.5)
render_asymptote("SecondOrderGradient.asy"; render=2)

Three points \$x,y,z\$ and the negative gradient of the second order absolute difference

If we now perform a gradient step, we obtain the three points

pn, rn, qn = exp.(Ref(M), [p, r, q], [-Xp, -Xr, -Xq])

as well we the new mid point

cn = mid_point(M, pn, qn)
geoPts_yncn = shortest_geodesic(M, rn, cn, T)
nothing #hide

and obtain the new situation

    points = [ [x,y,z], [c,c2,cn], [xn,yn,zn] ],
    curves = [ geoPts_yncn ] ,
    tVectors = [Tuple.([ [p, -Xp], [r, Xr], [q, Xq] ])],
    colors=Dict(:tvectors => [TolVibrantCyan],
        :points => [black, TolVibrantBlue, TolVibrantOrange],
        :curves => [TolVibrantTeal]
    dotSize = 3.5, lineWidth = 0.75, cameraPosition = (1.2,1.,.5)
render_asymptote("SecondOrderMin1.asy"; render=2)

A gradient Step

One can see, that this step slightly “overshoots”, i.e. $r$ is now even below $c$. and the cost function is still at

costTV2(M, (pn, rn, qn))

But we can also search for the best step size using linesearch_backtrack on the PowerManifold manifold $\mathcal N = \mathcal M^3 = (\mathbb S^2)^3$

x = [p, r, q]
N = PowerManifold(M, NestedPowerRepresentation(), 3)
s = linesearch_backtrack(
    x -> costTV2(M, Tuple(x)),
    [∇TV2(M, (p, r, q))...],  # transform from tuple to PowTVector
    1.0, # initial stepsize guess
    0.999, # decrease
    0.96,  #contract

and for the new points

pm, rm, qm = exp.(Ref(M), [p, r, q], s * [-Xp, -Xr, -Xq])
cm = mid_point(M, pm, qm)
geoPts_xmzm = shortest_geodesic(M, pm, qm, T)
nothing #hide

we obtain again with

    points = [ [x,y,z], [c,c2,cm], [xm,ym,zm] ],
    curves = [ geoPts_xmzm ] ,
    tVectors = [Tuple.( [-ξx, -ξy, -ξz], [x, y, z] )],
    colors=Dict(:tvectors => [TolVibrantCyan],
                :points => [black, TolVibrantBlue, TolVibrantOrange],
                :curves => [TolVibrantTeal]
    dotSize = 3.5, lineWidth = 0.75, cameraPosition = (1.2,1.,.5)

A gradient Step

Here, the cost function yields

costTV2(M, (pm, rm, qm))

which is nearly zero, as one can also see, since the new center $c$ and $r$ are quite close.