Lifting methods for manifold-valued variational problems

TL;DR

This paper reviews lifting methods that transform complex manifold-valued variational problems into convex problems in higher-dimensional spaces, enabling efficient solutions and extending techniques to general convex regularization.

Contribution

It introduces a refined finite element discretization framework for manifold-valued lifting methods and generalizes total variation regularization techniques.

Findings

Extended lifting methods to manifold-valued problems.

Supported general convex regularization in the lifted framework.

Facilitated global optimization of complex variational problems.

Abstract

Lifting methods allow to transform hard variational problems such as segmentation and optical flow estimation into convex problems in a suitable higher-dimensional space. The lifted models can then be efficiently solved to a global optimum, which allows to find approximate global minimizers of the original problem. Recently, these techniques have also been applied to problems with values in a manifold. We provide a review of such methods in a refined framework based on a finite element discretization of the range, which extends the concept of sublabel-accurate lifting to manifolds. We also generalize existing methods for total variation regularization to support general convex regularization.

Equations

$$\underset{u:\Omega \to \mathcal{M}}{\min }F(u)$$

$$\underset{u:\Omega \to \Gamma }{\min }F(u),$$

$$\underset{v:\Omega \to X}{\min }\tilde{F}(v),$$

$$F(u)={\int }_{\Omega }f(x,u(x),\nabla u(x))dx.$$

$$F(u)={\int }_{\Omega }\rho (x,u(x))dx+\lambda \mathrm{TV}(u),$$

$$\tilde{F}(v)={\int }_{\Omega }⟨\rho (x,\cdot ),v(x)⟩dx+\lambda \widetilde{\mathrm{TV}}(v),$$

$$\widetilde{\mathrm{TV}}(v):=\sup \left\{{\int }_{\Omega }⟨{\mathrm{div}}_{x}p(x,\cdot ),v(x)⟩dx:p:\Omega \times \mathcal{M}\to \mathbb{R},\Vert {\nabla }_{z}p{\Vert }_{\infty }\le 1\right\}.$$

$$\Vert {\nabla }_{z}p{\Vert }_{\infty }:=\sup \left\{\Vert {\nabla }_{z}p(x,z){\Vert }_{\sigma ,\infty }:(x,z)\in \Omega \times \mathcal{M}\right\},$$

$$Du(x)\xi \in T_{u(x)}\mathcal{M}\subset T_{u(x)}{\mathbb{R}}^N.$$

$$F(u):={\int }_{\Omega }f(x,u(x),Du(x))dx,$$

$$\tilde{F}(v):=\sup \left\{{\int }_{\Omega }⟨-{\mathrm{div}}_{x}p(x,\cdot )+q(x,\cdot ),v(x)⟩dx:({\nabla }_{z}p,q)\in \mathcal{K}\right\},$$

$$\mathcal{K}:=\left\{(P,q)\in C(\Sigma ;{\mathbb{R}}^{N,d}\times \mathbb{R}):f^{\ast }({\pi }_{1}y,{\pi }_{2}y,P(y))+q(y)\le 0\forall y\in \Sigma \right\},$$

$$f(x,z,\xi )=\rho (x,z)+\eta (P_z\xi )$$

$$\bigcup _{T\in {\mathcal{T}}_h}T={\mathcal{M}}_h.$$

$$S_h:=\{{\varphi }_h\in C^0({\mathcal{M}}_h):{\varphi }_h{|}_T\text{ is linear affine for each }T\in {\mathcal{T}}_h\}.$$

$${\mathbb{P}}_h({\mathcal{M}}_h)=\left\{{\mu }_h\in {\mathfrak{M}}_h({\mathcal{M}}_h):{\mu }_h\ge 0,\sum _{k=1}^L⟨{\mu }_h,{\chi }_k⟩=1\right\}.$$

$$\tilde{f}(x,z,\xi ):=\rho (x,\iota (z))+\eta (P_T\xi )$$

$$\forall x\in \Omega \forall z\in {\mathcal{M}}_h:{\eta }^{\ast }(P_T{\nabla }_{z}p(x,z))-\rho (x,\iota (z))+q(x,z)\le 0,$$

$${\eta }^{\ast }(P_T{\nabla }_{T}p(x))+⟨q_{T,1}(x),z⟩-\rho (x,\iota (z))\le -q_{T,2}(x),$$

$${\rho }_T^{\ast }(x,z):=\underset{z^{\prime }\in T}{\sup }⟨z,z^{\prime }⟩-\rho (x,\iota (z^{\prime })),$$

$${\eta }^{\ast }(P_T{\nabla }_{T}p(x))$$

$${\rho }_T^{\ast }(q_{T,1}(x))$$

$$a_T(x)+b_T(x)$$

$$\underset{v}{\inf }\underset{p,q}{\sup }$$

$${\eta }^{\ast }(P_T{\nabla }_{T}p(x))\le a_T(x),$$

$${\rho }_T^{\ast }(q_{T,1}(x))\le b_T(x),$$

$$a_T(x)+b_T(x)+q_{T,2}(x)=0.$$

$$\Vert P_T{\nabla }_{T}p(x){\Vert }_{\sigma ,\infty }\le \lambda \text{ and }q(x,z)\le \rho (x,\iota (z)).$$

$${\int }_{\Omega }\sum _{k=1}^L\rho (x,Z^k)v(x{)}^{k}dx.$$

$${\nabla }_{T}p=(p_1-p_2)\frac{Z_T^1-Z_T^2}{\Vert Z_T^1-Z_T^2{\Vert }_2^2}\in {\mathbb{R}}^2$$

$$P_T$$

$${\mu }_h=({\lambda }_1,\dots ,{\lambda }_L)\mapsto \underset{z\in \mathcal{M}}{\mathrm{argmin}}\sum _{k=1}^L{\lambda }_k{d}_{\mathcal{M}}(z,Z^k{)}^2$$

$$V_i^k$$

Full text

1 Introduction

Consider a variational image processing or general data analysis problem of the form

[TABLE]

with $\Omega\subset\mathbb{R}^{d}$ open and bounded. In this chapter, we will be concerned with problems where the image $u$ takes values in an $s$-dimensional manifold $\mathcal{M}$. Problems of this form are wide-spread in image processing and especially in the processing of manifold-valued images such as InSAR massonnet1998_vogt , EBSD bachmann2011_vogt , DTI basser1994_vogt , orientational/positional rosman2011_vogt data or images with values in non-flat color spaces such as hue-saturation-value (HSV) or chromaticity-brightness (CB) color spaces chan2001_vogt .

They come with an inherent non-convexity, as the space of images $u\colon\Omega\to\mathcal{M}$ is generally non-convex, with few exceptions, such as if $\mathcal{M}$ is a Euclidean space, or if $\mathcal{M}$ is a Hadamard manifold, if one allows for the more general notion of geodesic convexity bacak2014_vogt ; bacak2016_vogt . Except for these special cases, efficient and robust convex numerical optimization algorithms therefore cannot be applied and global optimization is generally out of reach.

The inherent non-convexity of the feasible set is not only an issue of representation. Even for seemingly simple problems, such as the problem of computing the Riemannian center of mass for a number of points on the unit circle, it can affect the energy in surprisingly intricate ways, creating multiple local minimizers and non-uniqueness (Fig. 1). The equivalent operation in Euclidean space, computing the weighted mean, is a simple convex (even linear) operation, with a unique, explicit solution.

The problem of non-convexity is not unique to our setting, but rather ubiquitous in a much broader context of image and signal processing: amongst others, image segmentation, 3D reconstruction, image matching, optical flow and image registration, superresolution, inpainting, edge-preserving image restoration with the Mumford-Shah and Potts model, machine learning, and many statistically or physically motivated models involve intrinsically non-convex feasible sets or energies. When applied to such non-convex problems, local optimization strategies often get stuck in local minimizers.

In convex relaxation approaches, an energy functional is approximated by a convex one whose global optimum can be found numerically and whose minimizers lie within a small neighborhood around the actual solution of the problem. A popular convex relaxation technique that applies to a wide range of problems from image and signal processing is functional lifting. With this technique, the feasible set is embedded into a higher-dimensional space where efficient convex approximations of the energy functional are easier available.

Overview and contribution. In the following sections, we will give a brief introduction to the concept of functional lifting and explore its generalization to manifold-valued problems. Our aim is to provide a survey-style introduction to the area, therefore we will provide references and numerical experiments on the way. In contrast to prior work, we will explain existing results in an updated finite element-based framework. Moreover, we propose extensions to handle general regularizers other than the total variation on manifolds, and to apply the “sublabel-accurate” methods to manifold-valued problems.

1.1 Functional lifting in Euclidean spaces

The problem of finding a function $u\colon\Omega\to\Gamma$ that assigns a label $u(x)\in\Gamma$ from a discrete range $\Gamma$ to each point $x$ in a continuous domain $\Omega\subset\mathbb{R}^{d}$, while minimizing an energy function $F(u)$, is commonly called a continuous multi-label (or multi-class labeling) problem in the image processing community pock2008_vogt ; lellmann2009_vogt . The name comes from the interpretation of this setting as the continuous counterpart to the fully discrete problem of assigning to each vertex of a graph one of finitely many labels $\gamma_{1},\dots,\gamma_{L}$ while minimizing a given cost function greig1989_vogt ; calinescu1998_vogt ; kleinberg2002_vogt ; ishikawa2003_vogt .

The prototypical application of multi-labeling techniques is multi-class image segmentation, where the task is to partition a given image into finitely many regions. In this case, the label set $\Gamma$ is discrete and each label represents one of the regions so that $u^{-1}(\gamma)\subset\Omega$ is the region that is assigned label $\gamma$.

In the fully discrete setting, one way of tackling first-order multi-label problems is to look for good linear programming relaxations calinescu1998_vogt ; kleinberg2002_vogt ; ishikawa2003_vogt . These approaches were subsequently translated to continuous domains $\Omega$ for the two-class chan2006_vogt , multi-class zach2008_vogt ; pock2008_vogt ; lellmann2011c_vogt ; bae2011_vogt , and vectorial goldluecke2010_vogt case, resulting in non-linear, but convex, relaxations. By honoring the continuous nature of $\Omega$, they reduce metrication errors and improve isotropy strekalovskiy2011_vogt ; strekalovskiy2012_vogt ; goldluecke2013_vogt ; strekalovskiy2015_vogt , see lellmann2013b_vogt for a discussion and more references.

The general strategy, which we will also follow for the manifold-valued case, is to replace the energy minimization problem

[TABLE]

by a problem

[TABLE]

where $X$ is some “nice” convex set of larger dimension than $\Gamma$ with the property that there is an embedding $i\colon\Gamma\hookrightarrow X$ and $F(u)\approx\tilde{F}(i\circ u)$ in some sense whenever $u\colon\Omega\to\Gamma$.

In general, the lifted functional $\tilde{F}$ is chosen in such a way that it exhibits favorable (numerical or qualitative) properties compared with the original functional $F$ while being sufficiently close to the original functional so that minimizers of $\tilde{F}$ can be expected to have some recoverable relationship with global minimizers of $F$. Usually, $\tilde{F}$ is chosen to be convex when $F$ is not, which will make the problem amenable for convex optimization algorithms and allows to find a global minimizer of the lifted problem.

While current lifting strategies generally avoid local minimizers of the original problem, they are still an approximation and they are generally not guaranteed to find the global minimizers of the original problem.

A central difficulty is that some simplifications have to be performed in the lifting process in order to make it computationally feasible, which may lose information about the original problem. As a result, global minimizers $v\colon\Omega\to X$ of the lifted problem need not be in the image of $\Gamma$ under the embedding $i\colon\Gamma\hookrightarrow X$ and therefore are not directly associated with a function in the original space.

The process of projecting a solution back to the original space of functions $u\colon\Omega\to\Gamma$ is a difficult problem and, unless $\Gamma$ is scalar pock2010_vogt , the projection cannot be expected to be a minimizer of the original functional (see the considerations in federer1974_vogt ; lavenant2017_vogt ; vogt2019_vogt ). These difficulties may be related to the fact that the original problems are NP-hard cremers2012_vogt . As in the discrete labeling setting kleinberg2002_vogt , so-called rounding strategies have been proposed in the continuous case lellmann2012_vogt ; lellmann2011phd_vogt that come with an a priori bound for the relative gap between the minimum of the original functional and the value attained at the projected version of a minimizer to the lifted functional. For the manifold-valued case considered here, we are not aware of a similar result yet.

In addition to the case of a discrete range $\Gamma$, relaxation methods have been derived for dealing with a continuous (non-discrete) range, most notably the scalar case $\Gamma\subseteq\mathbb{R}$ alberti2003_vogt ; pock2010_vogt . They typically consider first-order energies that depend pointwise on $u$ and $\nabla u$ only:

[TABLE]

The equivalent problem class in the fully discrete setting consists of the energies with only unary (depending on one vertex’s label) and pairwise (depending on two vertices’ labels) terms.

For the problem (4), applying a strategy as in (2)–(3) comes with a substantial increase in dimensions. These relaxation approaches therefore have been called functional lifting, starting from the paper pock2009_vogt where the (non-convex) Mumford-Shah functional for edge-preserving image regularization and segmentation is lifted to a space of functions $v\colon\Omega\times\Gamma\to[0,1]$, $\Gamma\subset\mathbb{R}$. The authors use the special “step function” lifting $X=\{v\colon\Gamma\to[0,1]\}$ and $i(z^{*})=v$ with $v(z)=1$ if $z\leq z^{*}$ and [math] otherwise, which is only available in the scalar case

In this case, the integrand $f\colon\Omega\times\Gamma\times\mathbb{R}^{s,d}\to\mathbb{R}$ in (4) is assumed to be convex in the third component and nonnegative. The less restrictive property of polyconvexity has been shown to be sufficient windheuser2016_vogt ; mollenhoff2019_vogt , so that also minimal surface problems fit into this framework. The continuous formulations can be demonstrated pock2009_vogt ; mollenhoff2019_vogt to have strong connections with the method of calibrations alberti2003_vogt and with the theory of currents giaquinta1998_vogt .

In this paper, we will consider the more general case of $\Gamma=\mathcal{M}$ having a manifold structure. We will also restrict ourselves to first-order models. Only very recently, attempts at generalizing the continuous lifting strategies to models with higher-order regularization have been made – for regularizers that depend on the Laplacian loewenhauser2018_vogt ; vogt2019_vogt in case of vectorial ranges $\Gamma\subset\mathbb{R}^{s}$ and for the total generalized variation ranftl2013_vogt ; strecke2019_vogt in case of a scalar range $\Gamma\subset\mathbb{R}$. However, in contrast to the first-order theory, the higher-order models, although empirically useful, are still considerably less mathematically validated. Furthermore, we mention that there are models where the image domain $\Omega$ is replaced by a shape (or manifold) delaunoy2009_vogt ; bernard2017_vogt , which is beyond the scope of this survey.

1.2 Manifold-valued functional lifting

In this chapter, we will be concerned with problems where $\Gamma$ has a manifold structure. The first step towards applying lifting methods to such problems was an application to the restoration of cyclic data strekalovskiy2011b_vogt ; cremers2012_vogt with $\Gamma=\mathbb{S}^{1}$, which was later lellmann2013_vogt generalized for the case of total variation regularization to data with values in more general manifolds. In lellmann2013_vogt , the functional lifting approach is applied to a first-order model with total variation regularizer,

[TABLE]

for $u\colon\Omega\to\mathcal{M}$, where $\Gamma=\mathcal{M}$ is an $s$-dimensional manifold and $\rho\colon\Omega\times\mathcal{M}\to\mathbb{R}$ is a pointwise data discrepancy. The lifted space is chosen to be $X=\mathbb{P}(\mathcal{M})$, the space of Borel probability measures over $\mathcal{M}$, with embedding $i\colon\mathcal{M}\hookrightarrow\mathbb{P}(\mathcal{M})$, where $i(z):=\delta_{z}$ is the Dirac point measure with unit mass concentrated at $z\in\mathcal{M}$ (see Fig. 2). The lifted functional is

[TABLE]

where $\langle g,\mu\rangle:=\int_{\mathcal{M}}g\,d\mu$ for $g\in C(\mathcal{M})$ and $\mu\in\mathbb{P}(\mathcal{M})$. Furthermore,

[TABLE]

The Lipschitz constraint $\|\nabla_{z}p\|_{\infty}\leq 1$, where

[TABLE]

and $\|\cdot\|_{\sigma,\infty}$ the spectral (operator) norm, can be explained by a functional analytic perspective vogt2018_vogt on this lifting strategy: The lifted total variation functional is the vectorial total variation semi-norm for functions over $\Omega$ with values in a certain Banach space of measures. The topological dual space of this space of measures is the space of Lipschitz continuous functions over $\mathcal{M}$. However, this interpretation does not generalize easily to other regularizers. We will instead base our model for general convex regularizers on the theory of currents as presented in mollenhoff2019_vogt .

Sublabel accuracy.

While the above model comes with a fully continuous description, a numerical implementation requires the discretization of $\Omega$ as well as the range $\Gamma$. This introduces two possible causes for errors: metrication errors and label bias.

Metrication errors are artifacts related to the graph or grid representation of the spatial image domain $\Omega$, finite difference operators, and the choice of metric thereof. They manifest mostly in unwanted anisotropy, missing rotational invariance, or blocky diagonals. They constitute a common difficulty with all variational problems and lifting approaches klodt2008_vogt .

In contrast, label bias means that the discretization favors solutions that assume values at the chosen “labels” (discretization points) $Z^{1},\dots,Z^{L}$ in the range $\Gamma$ (see Fig. 3 and 4). This is very desirable for discrete $\Gamma$, but in the context of manifolds, severely limits accuracy and forces a suitably fine discretization of the range.

In more recent so-called sublabel-accurate approaches for scalar and vectorial ranges $\Gamma$, more emphasis is put on the discretization zach2012_vogt ; mollenhoff2016_vogt ; laude2016_vogt to get rid of label bias in models with total variation regularization, which allows to greatly reduce the number of discretizations points for the range $\Gamma$. In a recent publication mollenhoff2017_vogt , the gain in sublabel accuracy is explained to be caused by an implicit application of first-order finite elements on $\Gamma$ as opposed to previous approaches that can be interpreted as using zero-order elements, which naturally introduces label-bias. An extension of the sublabel-accurate approaches to arbitrary convex regularizers using the theory of currents was recently proposed in mollenhoff2019_vogt .

Motivated by these recent advances, we propose to extend the methods from lellmann2013_vogt for manifold-valued images to arbitrary convex regularizers, making use of finite element techniques on manifolds dziuk2013_vogt . This reduces label bias and thus the amount of labels necessary in the discretization.

1.3 Further related work

The methods proposed in this work are applicable to variational problems with values in manifolds of dimension $s\leq 3$. The theoretical framework applies to manifolds of arbitrary dimension, but the numerical costs increase exponentially for dimensions $4$ and larger.

An alternative is to use local optimization methods on manifolds. A reference for the smooth case is absil2009_vogt . For non-smooth energies, methods such as the cyclic proximal point, Douglas-Rachford, ADMM and (sub-)gradient descent algorithm have been applied to first and second order TV and TGV as well as Mumford-Shah and Potts regularization approaches in weinmann2014_vogt ; weinmann2015_vogt ; baust2016_vogt ; bergmann2016_vogt ; bredies2018_vogt ; bergmann2018b_vogt . These methods are generally applicable to manifolds of any dimension whose (inverse) exponential mapping can be evaluated in reasonable time and quite efficient in finding a local miminum, but can get stuck in local extrema. Furthermore, the use of total variation regularization in these frameworks is currently limited to anisotropic formulations; Tikhonov regularization was proposed instead for isotropic regularization weinmann2014_vogt ; bergmann2018c_vogt . An overview of applications, variational models and local optimization methods is given in bergmann2018c_vogt .

Furthermore, we mention that, beyond variational models, there exist statistical fletcher2012_vogt , discrete graph-based bergmann2018_vogt , wavelet-based storath2018_vogt , PDE-based chefdhotel2004_vogt and patch-based laus2017_vogt models for the processing and regularization of manifold-valued signals.

2 Submanifolds of $\mathbb{R}^{N}$

We formulate our model for submanifolds of $\mathbb{R}^{N}$ which is no restriction by the Whitney embedding theorem (lee2013_vogt, , Thm. 6.15). For an $s$-dimensional submanifold of $\mathbb{R}^{N}$ and $\Omega\subset\mathbb{R}^{d}$ open and bounded, differentiable functions $u\colon\Omega\to\mathcal{M}$ are regarded as a subset of differentiable functions with values in $\mathbb{R}^{N}$. For those functions, a Jacobian $Du(x)\in\mathbb{R}^{N,d}$ in the Euclidean sense exists that can be identified with the push-forward of the tangent space $T_{x}\Omega$ to $T_{u(x)}\mathcal{M}$, i.e., for each $x\in\Omega$ and $\xi\in\mathbb{R}^{d}=T_{x}\Omega$, we have

[TABLE]

On the other hand, for differentiable maps $p\colon\mathcal{M}\to\mathbb{R}^{d}$, there exists an extension of $p$ to a neighborhood of $\mathcal{M}\subset\mathbb{R}^{N}$ that is constant in normal directions and we denote by $\nabla p(z)\in\mathbb{R}^{N,d}$ the Jacobian of this extension evaluated at $z\in\mathcal{M}$. Since the extension is assumed to be constant in normal directions, i.e., $\nabla p(z)\zeta=0$ whenever $\zeta\in N_{z}\mathcal{M}$ (the orthogonal complement of $T_{z}\mathcal{M}$ in $\mathbb{R}^{N}$), this definition is independent of the choice of extension.

2.1 Calculus of Variations on submanifolds

In this section, we generalize the total variation based approach in lellmann2013_vogt to less restrictive first-order variational problems by applying the ideas from functional lifting of vectorial problems mollenhoff2019_vogt to manifold-valued problems. Most derivations will be formal; we leave a rigorous choice of function spaces as well as an analysis of well-posedness for future work. We note that theoretical work is available for the scalar-valued case in alberti2003_vogt ; pock2010_vogt ; bouchitte2018_vogt and for the vectorial and for selected manifold-valued cases in giaquinta1998_vogt .

We consider variational models on functions $u\colon\Omega\to\mathcal{M}$,

[TABLE]

for which the integrand $f\colon\Omega\times\mathcal{M}\times\mathbb{R}^{N,d}\to\mathbb{R}$ is convex in the last component. Note that the dependence of $f$ on the full Jacobian of $u$ spares us dealing with the tangent bundle push-forward $T\Omega\to T\mathcal{M}$ in a coordinate-free way, thus facilitating discretization later on.

Formally, the lifting strategy for vectorial problems proposed in mollenhoff2019_vogt can be generalized to this setting by replacing the range $\Gamma$ with $\mathcal{M}$. As the lifted space, we consider the space of probability measures on the Borel $\sigma$-Algebra over $\mathcal{M}$, $X=\mathbb{P}(\mathcal{M})$, with embedding $i\colon\mathcal{M}\to\mathbb{P}(\mathcal{M})$, where $i(z)=\delta_{z}$ is the Dirac point mass concentrated at $z\in\mathcal{M}$. Furthermore, we write $\Sigma:=\Omega\times\mathcal{M}$ and, for $(x,z)=y\in\Sigma$, we define the coordinate projections $\pi_{1}y:=x$ and $\pi_{2}y:=z$. Then, for $v\colon\Omega\to\mathbb{P}(\mathcal{M})$, we define the lifted functional

[TABLE]

where $\langle g,\mu\rangle:=\int_{\mathcal{M}}g\,d\mu$ is the dual pairing between $g\in C(\mathcal{M})$ and $\mu\in\mathbb{P}(\mathcal{M})$ and

[TABLE]

where $f^{*}(x,z,\zeta):=\sup_{\xi}\langle\zeta,\xi\rangle-f(x,z,\xi)$ is the convex conjugate of $f$ with respect to the last variable.

In the following, the integrand $f\colon\Omega\times\mathcal{M}\times\mathbb{R}^{N,d}\to\mathbb{R}$ is assumed to decompose as

[TABLE]

into a pointwise data term $\rho\colon\Omega\times\mathcal{M}\to\mathbb{R}$ and a convex regularizer $\eta\colon\mathbb{R}^{s,d}\to\mathbb{R}$ that only depends on an $s$-dimensional representation of vectors in $T_{z}\mathcal{M}$ given by a surjective linear map $P_{z}\in\mathbb{R}^{s,N}$ with $\operatorname{ker}(P_{z})=N_{z}\mathcal{M}$.

This very general integrand covers most first-order models in the literature on manifold-valued imaging problems. It applies in particular to isotropic and anisotropic regularizers that depend on (matrix) norms of $Du(x)$ such as the Frobenius or spectral norm (or operator norm) where $P_{z}$ is taken to be an arbitrary orthogonal basis transformation. Since $z\mapsto P_{z}$ is not required to be continuous, it can also be applied to non-orientable manifolds such as the Moebius strip or the Klein bottle where no continuous orthogonal basis representation of the tangent bundle $T\mathcal{M}$ exists.

Regularizers of this particular form depend on the manifold through the choice of $P_{z}$ only. This is important because we approximate $\mathcal{M}$ in the course of our proposed discretization by a discrete (simplicial) manifold $\mathcal{M}_{h}$ and the tangent spaces $T_{z}\mathcal{M}$ are replaced by the linear spaces spanned by the simplicial faces of $\mathcal{M}_{h}$.

2.2 Finite elements on submanifolds

We translate the finite element approach for functional lifting proposed in mollenhoff2017_vogt to the manifold-valued setting by employing the notation from surface finite element methods dziuk2013_vogt .

The manifold $\mathcal{M}\subset\mathbb{R}^{N}$ is approximated by a triangulated topological manifold $\mathcal{M}_{h}\subset\mathbb{R}^{N}$ in the sense that there is a homeomorphism $\iota\colon\mathcal{M}_{h}\to\mathcal{M}$ (Fig. 5 and 6). By $\mathcal{T}_{h}$, we denote the set of simplices that make up $\mathcal{M}_{h}$:

[TABLE]

For $T,\tilde{T}\in\mathcal{T}_{h}$, either $T\cap\tilde{T}=\emptyset$ or $T\cap\tilde{T}$ is an ($s-k$)-dimensional face for $k\in\{1,\dots,s\}$. Each simplex $T\in\mathcal{T}_{h}$ spans an $s$-dimensional linear subspace of $\mathbb{R}^{N}$ and there is an orthogonal basis representation $P_{T}\in\mathbb{R}^{s,N}$ of vectors in $\mathbb{R}^{N}$ to that subspace. Furthermore, for later use, we enumerate the vertices of the triangulation as $Z^{1},\dots,Z^{L}\in\mathcal{M}\cap\mathcal{M}_{h}$.

For the numerics, we assume the first-order finite element space

[TABLE]

The functions in $S_{h}$ are piecewise differentiable on $\mathcal{M}_{h}$ and we define the surface gradient $\nabla_{T}\phi_{h}\in\mathbb{R}^{N,d}$ of $\phi_{h}\in S_{h}$ by the gradient of the linear affine extension of $\phi_{h}|_{T}$ to $\mathbb{R}^{N}$. If $L$ is the number of vertices in the triangulation of $\mathcal{M}_{h}$, then $S_{h}$ is a linear space of dimension $L$ with nodal basis $\chi_{1},\dots,\chi_{L}$ which is uniquely determined by the property $\chi_{k}(Z^{l})=1$ if $k=l$ and $\chi_{k}(Z^{l})=0$ otherwise (Fig. 7).

The dual space of $S_{h}$, which we denote by $\mathfrak{M}_{h}(\mathcal{M}_{h})$, is a space of signed measures. We identify $\mathfrak{M}_{h}(\mathcal{M}_{h})=\mathbb{R}^{L}$ via dual pairing with the nodal basis $\chi_{1},\dots,\chi_{L}$, i.e., to each $\mu_{h}\in\mathfrak{M}_{h}(\mathcal{M}_{h})$ we associate the vector $(\langle\mu_{h},\chi_{1}\rangle,\dots,\langle\mu_{h},\chi_{L}\rangle)$. We then replace the space $\mathbb{P}(\mathcal{M})$ of probability measures over $\mathcal{M}$ by the convex subset

[TABLE]

The energy functional is then translated to the discretized setting by redefining the integrand $f$ on $\mathcal{M}_{h}$ for any $x\in\Omega$, $z\in\mathcal{M}_{h}$ and $\xi\in\mathbb{R}^{N,d}$ as

[TABLE]

The epigraphical constraints in $\mathcal{K}$ translate to

[TABLE]

for functions $p\in S_{h}^{d}$ and $q\in S_{h}$. The constraints can be efficiently implemented on each $T\in\mathcal{T}_{h}$ where $\nabla_{z}p$ is constant and $q(x,z)=\langle q_{T,1}(x),z\rangle+q_{T,2}(x)$ is linear affine in $z$:

[TABLE]

for any $x\in\Omega$, $T\in\mathcal{T}_{h}$ and $z\in T$. Following the approach in mollenhoff2017_vogt , we define

[TABLE]

and introduce auxiliary variables $a_{T},b_{T}$ to split the epigraphical constraint (19) into two epigraphical and one linear constraint for $x\in\Omega$ and $T\in\mathcal{T}_{h}$:

[TABLE]

The resulting optimization problem is described by the following saddle point form over functions $v\colon\Omega\to\mathbb{P}_{h}(\mathcal{M}_{h})$, $p\in C^{1}(\Omega,S_{h}^{d+1})$ and $q\in C(\Omega,S_{h})$:

[TABLE]

Finally, for the fully discrete setting, the domain $\Omega$ is replaced by a Cartesian rectangular grid with finite differences operator $\nabla_{x}$ and Neumann boundary conditions.

2.3 Relation to lellmann2013_vogt

In lellmann2013_vogt , a similar functional lifting is proposed for the special case of total variation regularization and without the finite elements interpretation. More precisely, the regularizing term is chosen to be $\eta(\xi)=\lambda\|\xi\|_{\sigma,1}$ for $\xi\in\mathbb{R}^{s,d}$, where $\|\cdot\|_{\sigma,1}$ is the matrix nuclear norm, also known as Schatten-$1$-norm, which is given by the sum of singular values of a matrix. It is the dual to the matrix operator or spectral norm $\|\cdot\|_{\sigma,\infty}$. If we substitute this choice of $\eta$ into the discretization given above, the epigraphical constraint (18) translates to the two constraints

[TABLE]

The first one is a Lipschitz constraint just as in the model from lellmann2013_vogt , but two differences remain:

In lellmann2013_vogt , the lifted and discretized form of the data term reads

[TABLE]

This agrees with our setting if $z\mapsto\rho(x,\iota(z))$ is affine linear on each simplex $T\in\mathcal{T}_{h}$, as then $q(x,z)=\rho(x,\iota(z))$ maximizes the objective function for any $p$ and $v$. Hence, the model in lellmann2013_vogt doesn’t take into account any information about $\rho$ below the resolution of the triangulation. We improve this by implementing the epigraph constraints $\rho_{T}^{*}(q_{T,1}(x))\leq b_{T}(x)$ as proposed in laude2016_vogt using a convex approximation of $\rho_{T}$ (see Fig. 8). The approximation is implemented numerically with piecewise affine linear functions in a “sublabel-accurate” way, i.e., at a resolution below the resolution of the triangulation . 2. 2.

A very specific discretization of the gradients $\nabla_{T}p(x)$ is proposed in lellmann2013_vogt : To each simplex in the triangulation a mid-point $y_{T}\in\mathcal{M}$ is associated. The vertices $Z_{T}^{1},...,Z_{T}^{s+1}$ of the simplex are projected to the tangent space at $y_{T}$ as $v_{T}^{k}:=\log_{y_{T}}Z_{T}^{k}$. The gradient is then computed as the vector $g$ in the tangent space $T_{y_{T}}\mathcal{M}$ describing the affine linear map on $T_{y_{T}}\mathcal{M}$ that takes values $p(Z_{T}^{k})$ at the points $v_{T}^{k}$, $k=1,\dots,s+1$.

This procedure aims to make up for the error introduced by the simplicial discretization and amounts to a different choice of $P_{T}$ – a slight variant of our model. We did not observe any significant positive or negative effects from using either discretization; the difference between the minimizers is very small.

In the one-dimensional case, the two approaches differ only in a constant factor: Denote by $P_{T}\in\mathbb{R}^{s,N}$ the orthogonal basis representation of vectors in $\mathbb{R}^{N}$ in the subspace spanned by the simplex $T\in\mathcal{T}_{h}$ and denote by $\tilde{P}_{T}\in\mathbb{R}^{s,N}$ the alternative approach from lellmann2013_vogt . Now, consider a triangulation $\mathcal{T}_{h}$ of the circle $\mathbb{S}^{1}\subset\mathbb{R}^{2}$ and a one-dimensional simplex $T\in\mathcal{T}_{h}$. A finite element $p\in S_{h}$ that takes values $p_{1},p_{2}\in\mathbb{R}$ at the vertices $Z_{T}^{1},Z_{T}^{2}\in\mathbb{R}^{2}$ that span $T$ has the gradient

[TABLE]

and $P_{T},\tilde{P}_{T}\in\mathbb{R}^{1,2}$ are given by

[TABLE]

Hence $P_{T}=\alpha_{T}\tilde{P}_{T}$ for $\alpha_{T}=d_{\mathbb{S}^{1}}(Z_{T}^{1},Z_{T}^{2})/\|Z_{T}^{1}-Z_{T}^{2}\|_{2}$ the ratio between geodesic (angular) and Euclidean distance between the vertices. If the vertices are equally spaced on $\mathbb{S}^{1}$, this is a constant factor independent of $T$ that typically scales the discretized regularizer by a small constant factor. On higher-dimensional manifolds, more general linear transformations $P_{T}=A_{T}\tilde{P}_{T}$ come into play. For very irregular triangulations and coarse discretization, this may affect the minimizer; however, in our experiments the observed differences were negligible.

2.4 Full discretization and numerical implementation

A prime advantage of the lifting method when applied to manifold-valued problems is that it translates most parts of the problem into Euclidean space. This allows to apply established solution strategies for the non-manifold case, which rely on non-smooth convex optimization: After discretization, the convex-concave saddle-point form allows for a solution using the primal-dual hybrid gradient method chambolle2011_vogt ; chambolle2012_vogt with recent extensions goldstein2013_vogt . In this optimization framework, the epigraph constraints are realized by projections onto the epigraphs in each iteration step. For the regularizers to be discussed in this paper (TV, quadratic and Huber), we refer to the instructions given in pock2010_vogt . For the data term $\rho$, we follow the approach in laude2016_vogt : For each $x\in\Omega$, The data term $z\mapsto\rho(x,\iota(z))$ is sampled on a subgrid of $\mathcal{M}_{h}$ and approximated by a piecewise affine linear function. The quickhull algorithm can then be used to get the convex hull of this approximation. Projections onto the epigraph of $\rho^{*}_{T}$ are then projections onto convex polyhedra, which amounts to solving many low-dimensional quadratic programs; see laude2016_vogt for more details.

Following lellmann2013_vogt , the numerical solution $u\colon\Omega\to\mathbb{P}_{h}(\mathcal{M}_{h})$, taking values in the lifted space $\mathbb{P}_{h}(\mathcal{M}_{h})$, is projected back to a function $u\colon\Omega\to\mathcal{M}$, taking values in the original space $\mathcal{M}$, by mapping, for each $x\in\Omega$ separately, a probability measure $u(x)=(\lambda_{1},\dots,\lambda_{L})=\mu_{h}\in\mathbb{P}_{h}(\mathcal{M}_{h})$ to the following Riemannian center of mass on the original manifold $\mathcal{M}$:

[TABLE]

For $\mathcal{M}=\mathbb{R}^{s}$, this coincides with the usual weighted mean $\bar{z}=\sum_{k=1}^{L}\lambda_{k}Z_{k}$. However, on manifolds this minimization is known to be a non-convex problem with non-unique solutions (compare Fig. 1). Still, in practice the iterative method described in karcher1977_vogt yields reasonable results for all real-world data considered in this work: Starting from a point $z_{0}:=Z^{k}$ with maximum weight $\lambda_{k}$, we proceed for $i\geq 0$ by projecting the $Z^{k}$, $k=1,\dots,L$, to the tangent space at $z_{i}$ using the inverse exponential map, taking the linear weighted mean $v_{i}$ there and defining $z_{i+1}$ as the projection of $v_{i}$ to $\mathcal{M}$ via the exponential map:

[TABLE]

The method converges rapidly in practice. It has to be applied only once for each $x\in\Omega$ after solving the lifted problem, so that efficiency is non-critical.

3 Numerical Results

We apply our model to problems with quadratic data term $\rho(x,z):=d_{\mathcal{M}}^{2}(I(x),z)$ and Huber, total variation (TV) and Tikhonov (quadratic) regularization with parameter $\lambda>0$:

[TABLE]

where the Huber function $\phi_{\alpha}$ for $\alpha>0$ is defined by

[TABLE]

Note that previous lifting approaches for manifold-valued data were restricted to total variation regularization $\eta_{\text{TV}}$.

The methods were implemented in Python 3 with NumPy and PyCUDA, running on an Intel Core i7 4.00 GHz with an NVIDIA GeForce GTX 1080 Ti 12 GB and 16 GB RAM. The iteration was stopped as soon as the relative gap between primal and dual objective fell below $10^{-5}$. Approximate runtimes ranged between 5 and 45 minutes. The code is available from https://github.com/room-10/mfd-lifting.

3.1 One-dimensional denoising on a Klein bottle

Our model can be applied to both orientable and non-orientable manifolds. Figure 10 shows an application of our method to Tikhonov denoising of a synthetic one-dimensional signal $u\colon[0,1]\to\mathcal{M}$ on the two-dimensional Klein surface embedded in $\mathbb{R}^{3}$, a non-orientable closed surface that cannot be embedded into $\mathbb{R}^{3}$ without self-intersections. Our numerical implementation uses a triangulation with a very low count of $5\times 5$ vertices and 50 triangles. The resolution of the signal (250 one-dimensional data points) is far below the resolution of the triangulation and, still, our approach is able to restore a smooth curve.

3.2 Three-dimensional manifolds: $SO(3)$

Signals with rotational range $u\colon\Omega\to SO(3)$ occur in the description of crystal symmetries in EBSD (Electron Backscatter Diffraction Data) and in motion tracking. The rotation group $SO(3)$ is a three-dimensional manifold that can be identified with the three-dimensional unit-sphere $\mathbb{S}^{3}$ up to identification of antipodal points via the quaternion representation of 3D rotations. A triangulation of $\mathbb{S}^{3}$ is given by the vertices and simplicial faces of the hexacosichoron (600-cell), a regular polytope in $\mathbb{R}^{4}$ akin to the icosahedron in $\mathbb{R}^{3}$. As proposed in lellmann2013_vogt , we eliminate opposite points in the hexacosichoron and obtain a discretization of $SO(3)$ with $60$ vertices and $300$ tetrahedral faces.

Motivated by Bézier surface interpolation absil2016_vogt , we applied Tikhonov regularization to a synthetic inpainting (interpolation) problem with added noise (Fig. 11). In our variational formulation, we chose $\rho(x,z)=0$ for $x$ in the inpainting area and $\rho(x,z)=\delta_{\{z=I(x)\}}$ (a hard constraint to the input signal $I\colon\Omega\to SO(3)$) for $x$ in the known area.

Using the proposed sublabel-accurate handling of data terms, we obtain good results with only $60$ vertices, in contrast to lellmann2013_vogt , where the discretization is refined to $720$ vertices (Fig. 11).

3.3 Normals fields from digital elevation data

In digital elevation models (DEM), elevation information for earth science studies and mapping applications often includes surface normals which can be used to produce a shaded coloring of elevation maps. Normal fields $u\colon\Omega\to\mathbb{S}^{2}$ are defined on a rectangular image domain $\Omega\subset\mathbb{R}^{2}$; variational processing of the normal fields is therefore a manifold-valued problem on the two-dimensional sphere $\mathbb{S}^{2}\subset\mathbb{R}^{3}$.

Denoising using variational regularizers from manifold-valued image processing before computing the shading considerably improves visual quality (Fig. 12). For our framework, the sphere was discretized using 12 vertices and 20 triangles, chosen to form a regular icosahedron. The same dataset was used in lellmann2013_vogt , where the proposed lifting approach required 162 vertices – and solving a proportionally larger optimization problem – in order to produce comparable results.

We applied our approach with TV, Huber and Tikhonov regularization. Interestingly, many of the qualitative properties known from RGB and grayscale image processing appear to transfer to the manifold-valued case: TV enforces piecewise constant areas (flat hillsides), but preserves edges (mountain ridges). Tikhonov regularization gives overall very smooth results, but tends to lose edge information. With Huber regularization, edges (Mountain ridges) remain sharp while hillsides are smooth, and flattening is avoided (Fig. 12).

3.4 Denoising of high resolution InSAR data

While the resolution of the DEM dataset is quite limited ($40\times 40$ data points), an application to high resolution ($432\times 426$ data points) Interferometric Synthetic Aperture Radar (InSAR) denoising shows that our model is also applicable in a more demanding scenario (Fig. 13).

In InSAR imaging, information about terrain is obtained from satellite or aircraft by measuring the phase difference between the outgoing signal and the incoming reflected signal. This allows a very high relative precision, but no immediate absolute measurements, as all distances are only recovered modulo the wavelength. After normalization to $[0,2\pi)$, the phase data is correctly viewed as lying on the one-dimensional unit sphere $\mathbb{S}^{1}$. Therefore, handling the data before any phase unwrapping is performed requires a manifold-valued framework.

Again, denoising with TV, Huber, and Tikhonov regularizations demonstrates properties comparable to those known from scalar-valued image processing while all regularization approaches reduce noise substantially (Fig. 13).

4 Conclusion and Outlook

We provided an overview and framework for functional lifting techniques for the variational regularization of functions with values in arbitrary Riemannian manifolds. The framework is motivated from the theory of currents and continuous multi-label relaxations, but generalizes these from the context of scalar and vectorial ranges to geometrically more challenging manifold ranges.

Using this approach, it is possible to solve variational problems for manifold-valued images that consist of a possibly non-convex data term and an arbitrary, smooth or non-smooth, convex first-order regularizer, such as Tikhonov, total variation or Huber. A refined discretization based on manifold finite element methods achieves sublabel-accurate results, which allows to use coarser discretization of the range and reduces computational effort compared to previous lifting approaches on manifolds.

A primary limitation of functional lifting methods, which equally applies to manifold-valued models, is dimensionality: The numerical cost increases exponentially with the dimensionality of the manifold due to the required discretization of the range. Addressing this issue appears possible, but will require a significantly improved discretization strategy.

Acknowledgements.

The authors acknowledge support through DFG grant LE 4064/1-1 “Functional Lifting 2.0: Efficient Convexifications for Imaging and Vision” and NVIDIA Corporation.