# 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_field`

This example also illustrates the `PowerManifold`

manifold as well as `ArmijoLinesearch`

.

`using Manopt, Manifolds, Colors, PlutoUI`

We define some colors from Paul Tol

```
begin
black = RGBA{Float64}(colorant"#000000")
TolVibrantBlue = RGBA{Float64}(colorant"#0077BB")
TolVibrantOrange = RGBA{Float64}(colorant"#EE7733")
TolVibrantMagenta = RGBA{Float64}(colorant"#EE3377")
TolVibrantCyan = RGBA{Float64}(colorant"#33BBEE")
TolVibrantTeal = RGBA{Float64}(colorant"#009988")
end;
```

```
begin
T = [0:0.1:1.0...]
#render asy yes/no. If not, images included w/ markdown are assumed to be prerendered
render_asy = false
localpath = join(splitpath(@__FILE__)[1:(end - 1)], "/") # remove file to get files folder
image_prefix = localpath * "/second_order_difference"
@info image_prefix
end;
```

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

```
begin
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)
end;
```

Now the second order absolute difference can be stated as [BacakBergmannSteidlWeinmann2016]

$$d_2(p_1,p:2,p_3) := \min_{c ∈ \mathcal C_{p_1,p_3}} d_{\mathcal M}(c,p_2),\qquad p_1,p_2,p_3∈\mathcal M,$$

where $\mathcal C_{p,q}$ is the set of all mid points $g(\frac{1}{2};p,q)$, between `p`

and `q`

, i.e. where $g$ is a (not necessarily minimizing) geodesic connecting both.

For illustration we further define the point opposite of the mid point `c`

defined above

`c2 = -c;`

To illustrate the second order difference let‘s look at the geodesic connecting $r$ and the mid point $c$

`geoPts_rc = shortest_geodesic(M, r, c, T);`

```
render_asy && begin
asymptote_export_S2_signals(
image_prefix * "/SecondOrderData.asy";
curves=[geoPts_rc],
points=[[p, r, q], [c, c2]],
colors=Dict(:curves => [TolVibrantTeal], :points => [black, TolVibrantBlue]),
dot_size=3.5,
line_width=0.75,
camera_position=(1.2, 1.0, 0.5), #src
)
render_asymptote(image_prefix * "/SecondOrderData.asy"; render=2)
end;
```

`PlutoUI.LocalResource(image_prefix * "/SecondOrderData.png")`

**Figure.** *The three poins $p, q$, and $r$ (neart north pole, all black), which is connected by a geodesic (teal) to the mid point of the first two, $c$, (blue). The length of this geodesic is the cost of the second order total variation.*

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 the second order total variation cost, see `costTV2`

, yields.

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

`1.413716694115407`

But also its gradient can be easily computed since it is just a distance with respect to $r$ 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. The gradient is also directly implemented, see `grad_TV2`

. 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]
```

```
render_asy && begin
asymptote_export_S2_signals(
image_prefix * "/SecondOrderGradient.asy";
points=[[p, r, q], [c, c2]],
tangent_vectors=[Tuple.([[p, -Xp], [r, -Xr], [q, -Xq]])],
colors=Dict(:tvectors => [TolVibrantCyan], :points => [black, TolVibrantBlue]),
dot_size=3.5,
line_width=0.75,
camera_position=(1.2, 1.0, 0.5),
)
render_asymptote(image_prefix * "/SecondOrderGradient.asy"; render=2)
end;
```

`PlutoUI.LocalResource(image_prefix * "/SecondOrderGradient.png")`

**Figure.** *The negative gradient of the second order difference cost indicates the movement of the three points in order to reduce their cost.*

If we now perform a gradient step with constant step size 1, 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)`

```
3-element Vector{Float64}:
0.4508003528726723
0.4508003528726723
0.7704272085666162
```

Let‘s also again consider the geodesic connecting the new point $r_n$ and the new mid point $c_n$ as well as the gradient

```
begin
geoPts_rncn = shortest_geodesic(M, rn, cn, T)
(Xpn, Xrn, Xqn) = grad_TV2(M, (pn, rn, qn))
end;
```

The new configuration of the three points looks as follows

```
render_asy && begin
asymptote_export_S2_signals(
image_prefix * "/SecondOrderMin1.asy";
points=[[p, r, q], [c, c2, cn], [pn, rn, qn]],
curves=[geoPts_rncn],
tangent_vectors=[
Tuple.([[p, -Xp], [r, -Xr], [q, -Xq]]),
Tuple.([[pn, -Xpn], [rn, -Xrn], [qn, -Xqn]]),
],
colors=Dict(
:tvectors => [TolVibrantCyan, TolVibrantMagenta],
:points => [black, TolVibrantBlue, TolVibrantOrange],
:curves => [TolVibrantTeal],
),
dot_size=3.5,
line_width=0.75,
camera_position=(1.2, 1.0, 0.5),
)
render_asymptote(image_prefix * "/SecondOrderMin1.asy"; render=2)
end;
```

`PlutoUI.LocalResource(image_prefix * "/SecondOrderMin1.png")`

**Figure.** *The new situation of $p_n, q_n$, and $r_n$ (orange) and the miod point of the first two, $c_n$ (blue), which is again connected by a geodesic to $r_n$. Note that this geodesic is shorter, but also that $c$ and $r$ switched places*. The new gradient (magenta) is only slightly reduced in magnitude.*

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))`

`0.46579428818288593`

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$

```
begin
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.957, 0.999)
end
```

`0.7362738857656108`

and we obtain the new points

```
begin
pm, rm, qm = exp.(Ref(M), [p, r, q], s * [-Xp, -Xr, -Xq])
cm = mid_point(M, pm, qm)
geoPts_pmqm = shortest_geodesic(M, pm, qm, T)
end;
```

we obtain

```
render_asy && begin
asymptote_export_S2_signals(
image_prefix * "/SecondOrderMin2.asy";
points=[[p, r, q], [c, c2, cm], [pm, rm, qm]],
curves=[geoPts_pmqm],
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.0, 0.5),
)
render_asymptote(image_prefix * "/SecondOrderMin2.asy"; render=2)
end;
```

`PlutoUI.LocalResource(image_prefix * "/SecondOrderMin2.png")`

**Figure.** *For the best step size found by line search, $r_m$ and $c_m$ nearly agree, i.e. $r_m$ lies on the geodesic between $p_m$ and $q_m$ as the geodesic drawn here indicates.*

Here, the cost function yields

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

`0.003907643132786784`

which is much closer to zero, as one can also see, since the new center $c_m$ and $r_m$ are quite close.

## Literature

BacakBergmannSteidlWeinmann2016

Bačák, M., Bergmann, R., Steidl, G. and Weinmann, A.:

A second order nonsmooth variational model for restoring manifold-valued images, SIAM Journal on Scientific Computations, Volume 38, Number 1, pp. A567–597, doi: 10.1137/15M101988X, arXiv: 1506.02409