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

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]
2-element Vector{Vector{Float64}}:
 [0.7071067811865475, 0.7071067811865475, 0.0]
 [0.11061587104123713, 0.11061587104123713, 0.9876883405951378]

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
3-element Vector{Float64}:

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)

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]),
    dot_size = 3.5, line_width = 0.75, camera_position = (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 grad_TV2 for details. we obtain

(Xp, Xr, Xq) = grad_TV2(M, (p, r, q))
3-element Vector{Vector{Float64}}:
 [-0.0, -4.9676995583751974e-18, -0.7071067811865475]
 [-0.6984011233337103, -0.6984011233337102, 0.15643446504023084]
 [4.9676995583751974e-18, 0.0, -0.7071067811865475]

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]),
   dot_size = 3.5, line_width = 0.75, camera_position = (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])
3-element Vector{Vector{Float64}}:
 [0.7602445970756302, 4.563951614149274e-18, 0.6496369390800624]
 [0.6474502912317517, 0.6474502912317516, 0.4020152245473301]
 [-4.563951614149274e-18, 0.7602445970756302, 0.6496369390800624]

as well we the new mid point

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

and obtain the new situation

    points = [ [x,y,z], [c,c2,cn], [xn,yn,zn] ],
    curves = [ geoPts_yncn ] ,
    tangent_vectors = [Tuple.([ [p, -Xp], [r, -Xr], [q, -Xq] ])],
    colors=Dict(:tvectors => [TolVibrantCyan],
        :points => [black, TolVibrantBlue, TolVibrantOrange],
        :curves => [TolVibrantTeal]
    dot_size = 3.5, line_width = 0.75, camera_position = (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(N, x -> costTV2(N, x), x, grad_TV2(N, x), 1.0, 0.96, 0.999)

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)

we obtain again with

    points = [ [x,y,z], [c,c2,cm], [xm,ym,zm] ],
    curves = [ geoPts_xmzm ] ,
    tangent_vectors = [Tuple.( [-ξx, -ξy, -ξz], [x, y, z] )],
    colors=Dict(:tvectors => [TolVibrantCyan],
                :points => [black, TolVibrantBlue, TolVibrantOrange],
                :curves => [TolVibrantTeal]
    dot_size = 3.5, line_width = 0.75, camera_position = (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.