Norwegian University of Science and Technology

Fenchel Duality Theory anda Primal-Dual Algorithm
on Riemannian Manifolds
Ronny Bergmann
joint work with

R. Herzog, M. Silva Louzeiro, D. Tenbrinck, J. Vidal-Núñez.
40nt International Conference on Machine Learning
Workshop on Duality Principles for Modern Machine Learning
Honolulu, Hawaii,

July 29, 2023
travel support by RIKEN

Manifold-valued Signal & Image Processing
Tasks in image processing are often phrased as an optimisation problem.
Here. The pixel take values on a manifold▶ phase-valued data (S1 )▶ wind-fields, GPS (S2 )▶ DT-MRI (P(3))▶ EBSD, (grain) orientations (SO(n))
Artificial noisy phase-valued data.

Tasks. Denoising, Inpainting, labeling (classification), deblurring,...
2

Manifold-valued Signal & Image Processing
Tasks in image processing are often phrased as an optimisation problem.
Here. The pixel take values on a manifold▶ phase-valued data (S1 )▶ wind-fields, GPS (S2 )▶ DT-MRI (P(3))▶ EBSD, (grain) orientations (SO(n))
InSAR-Data of Mt. Vesuvius.

[Rocca, Prati, and Guarnieri 1997]

Tasks. Denoising, Inpainting, labeling (classification), deblurring,...
2

Manifold-valued Signal & Image Processing
Tasks in image processing are often phrased as an optimisation problem.
Here. The pixel take values on a manifold▶ phase-valued data (S1 )▶ wind-fields, GPS (S2 )▶ DT-MRI (P(3))▶ EBSD, (grain) orientations (SO(n))
Artificial noisy data on the sphere S2 .

Tasks. Denoising, Inpainting, labeling (classification), deblurring,...
2

Manifold-valued Signal & Image Processing
Tasks in image processing are often phrased as an optimisation problem.
Here. The pixel take values on a manifold▶ phase-valued data (S1 )▶ wind-fields, GPS (S2 )▶ DT-MRI (P(3))▶ EBSD, (grain) orientations (SO(n))
Artificial diffusion data,
each pixel is a symmetric positive matrix.

Tasks. Denoising, Inpainting, labeling (classification), deblurring,...
2

Manifold-valued Signal & Image Processing
Tasks in image processing are often phrased as an optimisation problem.
Here. The pixel take values on a manifold▶ phase-valued data (S1 )▶ wind-fields, GPS (S2 )▶ DT-MRI (P(3))▶ EBSD, (grain) orientations (SO(n))
DT-MRI of the human brain.
Camino Profject: cmic.cs.ucl.ac.uk/camino

Tasks. Denoising, Inpainting, labeling (classification), deblurring,...
2

Manifold-valued Signal & Image Processing
Tasks in image processing are often phrased as an optimisation problem.
Here. The pixel take values on a manifold▶ phase-valued data (S1 )▶ wind-fields, GPS (S2 )▶ DT-MRI (P(3))▶ EBSD, (grain) orientations (SO(n))
Grain orientations in EBSD data.
MTEX toolbox: mtex-toolbox.github.ioTasks. Denoising, Inpainting, labeling (classification), deblurring,...
2

A Riemannian Manifold M

A d-dimensional Riemannian manifold can be informally defined as a set M
covered with a ‘suitable’ collection of charts, that identify subsets of M with
open subsets of Rd and a continuously varying inner product on the tangent[Absil, Mahony, and Sepulchre 2008]
spaces.
Tp M

Xlogp pexpp
logp

Notation.
▶ Logarithmic map logp q = γ̇(0; p, q)▶ Exponential map expp X = γp,X (1)▶ Geodesic γ(·; p , q)

Pq←p YYpγ(·; p , q)

▶ Tangent space Tp M
▶ inner product (· , ·)p3

▶ parallel transport Pq←p XMq

The Model
We consider a minimization problemarg min F(p) + G(Λ(p))
p∈C▶ M, N are (high-dimensional) Riemannian Manifolds▶ F : M → R nonsmooth, (locally, geodesically) convex▶ G : N → R nonsmooth, (locally) convex▶ Λ : M → N nonlinear▶ C ⊂ M strongly geodesically convex.
In image processing.
choose a model, such that finding a minimizer yields the reconstruction
4

Splitting Methods & Algorithms
On a Riemannian manifold M we have
▶ Cyclic Proximal Point Algorithm (CPPA)[Bačák 2014]▶ (parallel) Douglas–Rachford Algorithm (PDRA)
On Rn PDRA is known to be equivalent to[RB, Persch, and Steidl 2016][Setzer 2011; O’Connor and Vandenberghe 2018]▶ Primal-Dual Hybrid Gradient Algorithm (PDHGA)[Esser, Zhang, and Chan 2010]▶ Chambolle-Pock Algorithm (CPA)[Chambolle and Pock 2011; Pock, Cremers, Bischof, and Chambolle 2009]

But on a Riemannian manifold M: đ no duality theory!
Goals of this talk.
Formulate Duality on a Manifold
Derive a Riemannian Chambolle–Pock Algorithm (RCPA)
5The Euclidean Fenchel Conjugate
Let f : Rn → R be proper and convex.
We define the Fenchel conjugate f ∗ : Rn → R of f by(∗f (ξ) := sup hξ, xi − f(x) = supx∈Rnx∈Rnξ−1)T (xf(x)▶ interpretation: maximize the distance of ξT x to f
⇒ extremum seeking problem on the epigraph
The Fenchel biconjugate readsf ∗ ∗ (x) = (f ∗ ) ∗ (x) = sup hξ , xi − f ∗ (ξ).
ξ∈Rn6)

Illustration of the Fenchel Conjugate
The Fenchel conjugate f∗The function f∗
1 f (ξ)4 f(x)f∗(ξ)

30.5

2

ξ−41x−1−f∗(ξ)−1712−22−0.5−1The Riemannian m−Fenchel Conjugate[RB, Herzog, Silva Louzeiro, Tenbrinck, and Vidal-Núñez 2021]
alternative approaches: [Ahmadi Kakavandi and Amini 2010; Silva Louzeiro, RB, and Herzog 2022]

Idea: Introduce a point on M to “act as” 0.
Let m ∈ C ⊂ M be given and F : C → R.
The m-Fenchel conjugate F∗m : Tm∗M → R is defined by{}F∗m (ξm ) := sup hξm , Xi − F(expm X) ,X∈LC,mwhere LC,m := {X ∈ Tm M | q = expm X ∈ C and kXkp = d(q, p)}.
Let m′ ∈ C. The mm′ -Fenchel-biconjugate F∗∗
mm′ : C → R is given by}{hξm′ , logm′ pi − F∗m (Pm←m′ ξm′ ) .
supF∗∗mm′ (p) =ξm′ ∈Tm∗′M

usually we only use the case m = m′ .
8

Properties of the m-Fenchel Conjugate[RB, Herzog, Silva Louzeiro, Tenbrinck, and Vidal-Núñez 2021]▶ F∗m is convex on Tm∗M▶ F(p) ≤ G(p) for all p ∈ C ⇒ F∗m (ξm ) ≥ G∗m (ξm ) for all ξm ∈ Tm∗M
▶ Fenchel-Moreau theorem: F ◦ expm convex (on Tm M), proper, lsc,then F∗∗
mm = F on C.
▶ Fenchel-Young inequality: For a proper, convex function F ◦ expm
ξp ∈ ∂M F(p) ⇔ F(p) + F∗m (Pm←p ξp ) = hPm←p ξp , logm pi.
▶ For a proper, convex, lsc function F ◦ expm
ξp ∈ ∂M F(p) ⇔ logm p ∈ ∂F∗m (Pm←p ξp ).

9Proximal Map

For f : M → R and λ > 0 we define the Proximal Map as[Moreau 1965; Rockafellar 1970; Ferreira and Oliveira 2002]

proxλf (p) := arg min dM (u, p)2 + λf(u).
u∈M! For a minimizer u∗ of f we have proxλf (u∗ ) = u∗ .
▶ For f proper, convex, lsc:
▶ the proximal map is unique.
▶ Proximal-Point-Algorithm:
pk = proxλf (pk−1 ) converges to arg min f10

The Chambolle-Pock Algorithm[Chambolle and Pock 2011]

From the pair of primal-dual problems
min f(x) + g(Kx),x∈RnK linear,max − f ∗ (−K∗ ξ) − g∗ (ξ)ξ∈Rm

we obtain for f, g proper convex, lsc theˆ as ,optimality conditions (OC) for a solution (x̂, ξ)Chambolle–Pock Algorithm. with σ > 0, τ > 0, θ ∈ R reads∂f3 −K∗ ξˆˆ Kx̂
∂g∗ (ξ)3ξ¯(k+1) = ξ (k+1) + θ(ξ (k+1) − ξ (k) )11

The Exact Riemannian Chambolle–Pock Algorithm[RB, Herzog, Silva Louzeiro, Tenbrinck, and Vidal-Núñez 2021; Chambolle and Pock 2011]

Assume. f(p) = F(p) + G(Λ(p)), with Λ : M → N .
Input: m, p(0) ∈ C ⊂ M, n = Λ(m), ξn ∈ Tn∗N , and parameters σ, τ, θ > 0
1: k ← 0
2: p̄(0) ← p(0)3: while not converged do( (k)())
(k+1)
4:
ξn
← proxτ G∗n ξn + τ logn Λ(p̄(k) ) ♭(())((k+1) )♯(k+1)∗
5:p← proxσF expp(k) Pp(k) ←m − σDΛ(m) [ξn]()6:p̄(k+1) ← expp(k+1) −θ logp(k+1) p(k)7:k←k+1
8: end while
Output: p(k)(0)

12Generalizations & Variants of the RCPAClassically[Chambolle and Pock 2011]▶ change σ = σk , τ = τk , θ = θk during the iterations▶ introduce an acceleration γ▶ relax dual ξ¯ instead of primal p̄ (switches lines 4 and 5)
Furthermore we[RB, Herzog, Silva Louzeiro, Tenbrinck, and Vidal-Núñez 2021]▶ introduce the lRCPA: linearize Λ, i. e., adopt the Euclidean case from

[Valkonen 2014](k)logn Λ(p̄→)(k)

Pn←Λ(m) DΛ(m)[logm p̄]▶ choose n 6= Λ(m) introduces a parallel transportDΛ(m)∗ [ξn(k+1)]→DΛ(m)∗ [PΛ(m)←n ξn▶ change m = m(k) , n = n(k) during the iterations
13(k+1)]

ManifoldsBase.jl & Manifolds.jl
ManifoldsBase.jl is an interface
for Riemannian manifolds M▶ inner(M, p, X, Y) (X , Y)p▶ exp(M, p, X) and log(M, p, q),▶ more general:
retract(M, p, X, m),where m is a retraction method▶ embeddings as decoratorѸmutating variants, e. g.
exp!(M, q, p, X)

works in place of q

Manifolds.jl is a Library of manifolds▶ Circle, (unit) Sphere & Torus▶ Fixed Rank Matrices▶ (Symplectic) Stiefel & Grassmann▶ Hyperbolic space & Rotations▶ Symmetric positive definite matrices▶ ...and many more
as well as generically▶ power & product manifold▶ tangent & vector bundles▶ Lie groups, connections, metrics,…

14

Ń juliamanifolds.github.io/ManifoldsBase.jl/
Ń juliamanifolds.github.io/Manifolds.jl/[Axen, Baran, RB, and Rzecki 2021]

Manopt.jl: Optimisation on Manifolds in Julia
Goal. Optimisation algorithms on Riemannian manifolds, based onManifoldsBase.jl ⇒ works with any manifold from Manifolds.jl.
Features.
▶ generic algorithm framework:
With Problem p and a SolverState s▶ initialize_solver!(p, s)▶ step_solver!(p, s, i): ith step

ˣ run algorithm: call solve(p, s)▶ generic debug and recording▶ step sizes and stopping criteria.
Manopt Family.

15[RB 2022]
manoptjl.org[Boumal,Mishra,Absil,andSepulchre
2014]
manopt.org
pymanopt.org [Townsend, Koep, and Weichwald 2016]

Algoirthms.
▶ Nelder-Mead, Particle Swarm▶ Subgradient Method▶ Gradient Descent
CG, Stochastic, Momentum, ...

▶ Quasi-Newton
BFGS, DFP, Broyden, SR1, ...

▶ Trust Regions▶ Chambolle-Pock▶ Douglas-Rachford, CPPA
▶ ALM, EPM, Frank-Wolfe,...
▶ Difference of Convex
DCA, DCPPAThe ℓ2 -TV Model[Rudin, Osher, and Fatemi 1992; Lellmann, Strekalovskiy, Koetter, and Cremers 2013; Weinmann, Demaret, and Storath 2014]

For a manifold-valued image f ∈ M, M = N d1 ,d2 , we computearg minp∈M

1F(p) + G(Λ(p)),αα > 0,with▶ data term F(p) = 12 d2M (p, f)▶ “forward differences” Λ : M → (TM)d1 −1, d2 −1, 2 ,()p 7→ Λ(p) = (logpi pi+e1 , logpi pi+e2 )i∈{1,...,d1 −1}×{1,...,d2 −1}▶ prior G(X) = kXkg,q,1 similar to a collaborative TV[Duran, Moeller, Sbert, and Cremers 2016]⇒ proxλG∗n given in closed form for q = 1 (anisotropic TV) and q = 2(isotropic TV).
16

Numerical Example for a P(3)-valued Image

P(3)-valued data.
anisotropic TV, α = 6.
▶ in each pixel we have a symmetric positive definite matrix▶ Applications: denoising/inpainting e.g. of DT-MRI data
17

Numerical Example for a P(3)-valued Image

P(3)-valued data.

anisotropic TV, α = 6.

Approach. CPPA as benchmark
parameters
iterations
runtime17[Bačák 2014; RB, Herzog, Silva Louzeiro, Tenbrinck, and Vidal-Núñez 2021]

CPPA

PDRA

lRCPA

λk = 4k

λ = 0.58
β = 0.93

σ = τ = 0.4
γ = 0.2, m = I

4000
1235 s.

122
380 s.

113
96.1 s.

Numerical Example for a P(3)-valued Image
38.74
CPPA
PDRA
lRCPA

Cost

80
60
40
1

10Approach. CPPA as benchmark
parameters
iterations
runtime
17100

Iterations

1,000[Bačák 2014; RB, Herzog, Silva Louzeiro, Tenbrinck, and Vidal-Núñez 2021]

CPPAPDRAlRCPAλk = 4k

λ = 0.58
β = 0.93

σ = τ = 0.4
γ = 0.2, m = I

4000
1235 s.

122
380 s.

113
96.1 s.

SummarySummary.
▶ We introduced a duality framework on manifolds▶ we introduced a Riemannian Chambolle–Pock algorithm▶ We saw a Software framework for Optimisation algorithms on manifolds▶ Numerical examples illustrates its performanceˣ Another model works with both functions being geodesically convex[Silva Louzeiro, RB, and Herzog 2022]

ľ Outlook.

▶ Explore further areas where Duality can be used in non-Euclidean spaces▶ Explore further connections between Duality-based algorithms▶ look into further applications
18

Selected References
Õ Axen, S. D., M. Baran, RB, and K. Rzecki (2021). Manifolds.jl: An Extensible Julia Frameworkfor Data Analysis on Manifolds. arXiv: 2106.08777.

Õ Bačák, M. (2014). “Computing medians and means in Hadamard spaces”. In: SIAM Journal on
Optimization 24.3, pp. 1542–1566. doi: 10.1137/140953393.

Õ RB (2022). “Manopt.jl: Optimization on Manifolds in Julia”. In: Journal of Open Source Software
7.70, p. 3866. doi: 10.21105/joss.03866.

Õ RB, R. Herzog, M. Silva Louzeiro, D. Tenbrinck, and J. Vidal-Núñez (Jan. 2021). “Fenchel
duality theory and a primal-dual algorithm on Riemannian manifolds”. In: Foundations ofComputational Mathematics. doi: 10.1007/s10208-020-09486-5. arXiv: 1908.02022.

Õ RB, J. Persch, and G. Steidl (2016). “A parallel Douglas Rachford algorithm for minimizing
ROF-like functionals on images with values in symmetric Hadamard manifolds”. In: SIAM Journal
on Imaging Sciences 9.4, pp. 901–937. doi: 10.1137/15M1052858.

Õ Chambolle, A. and T. Pock (2011). “A first-order primal-dual algorithm for convex problems withapplications to imaging”. In: Journal of Mathematical Imaging and Vision 40.1, pp. 120–145. doi:
10.1007/s10851-010-0251-1.

Õ Silva Louzeiro, M., RB, and R. Herzog (June 2022). “Fenchel Duality and a Separation Theorem
on Hadamard Manifolds”. In: SIAM Journal on Optimization 32.2, pp. 854–873. issn: 1052-6234,1095-7189. doi: 10.1137/21M1400699. arXiv: 2102.11155.

19

Ǥ ronnybergmann.net/talks/2023-ICML-Riemannian-Duality.pdf