Illustration of the Gradient of a Second Order Difference
This example explains how to compute the gradient of the second order difference midpoint model using adjoint_Jacobi_field
.
This example also illustrates the PowerManifold
manifold as well as ArmijoLinesearch
.
Setup
If you open this notebook in Pluto locally it switches between two modes. If the tutorial is within the Manopt.jl
repository, this notebook tries to use the local package in development mode. Otherwise, the file uses the Pluto pacakge management version. In this case, the includsion of images might be broken. unless you create a subfolder optimize
and activate asy
-rendering.
Since the loading is a little complicated, we show, which versions of packages were installed in the following.
with_terminal() do
Pkg.status()
end
�[32m�[1mStatus�[22m�[39m `/private/var/folders/_v/wg192lpd3mb1lp55zz7drpcw0000gn/T/jl_NnARuh/Project.toml` �[90m [5ae59095] �[39mColors v0.12.10 �[90m [1cead3c2] �[39mManifolds v0.8.42 �[90m [0fc0a36d] �[39mManopt v0.4.0 `~/Repositories/Julia/Manopt.jl` �[90m [7f904dfe] �[39mPlutoUI v0.7.49 �[90m [44cfe95a] �[39mPkg v1.8.0 �[90m [9a3f8284] �[39mRandom
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"
_in_package && @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 midpoints $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 midpoint c
defined above
c2 = -c;
To illustrate the second order difference let's look at the geodesic connecting $r$ and the midpoint $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;
Figure.The three (black) points $p, q$, and $r$ (near the north pole), which is connected by a geodesic (teal) to the midpoint 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_shortest_geodesic_startpoint
and adjoint_differential_shortest_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;
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 midpoint
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 midpoint $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 midpoint 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