Lifting Vectorial Variational Problems: A Natural Formulation based on Geometric Measure Theory and Discrete Exterior Calculus

TL;DR

This paper introduces a novel convex formulation for vectorial variational problems in imaging by lifting them to the space of currents and discretizing with Whitney forms, enabling more effective shape optimization.

Contribution

It presents a new convex relaxation approach for vector-valued variational problems using geometric measure theory and discrete exterior calculus, extending multilabeling methods.

Findings

Convex relaxation via currents improves problem tractability.

Discretization with Whitney forms generalizes multilabeling approaches.

The method facilitates shape optimization in imaging tasks.

Abstract

Numerous tasks in imaging and vision can be formulated as variational problems over vector-valued maps. We approach the relaxation and convexification of such vectorial variational problems via a lifting to the space of currents. To that end, we recall that functionals with polyconvex Lagrangians can be reparametrized as convex one-homogeneous functionals on the graph of the function. This leads to an equivalent shape optimization problem over oriented surfaces in the product space of domain and codomain. A convex formulation is then obtained by relaxing the search space from oriented surfaces to more general currents. We propose a discretization of the resulting infinite-dimensional optimization problem using Whitney forms, which also generalizes recent "sublabel-accurate" multilabeling approaches.

Figures

Equations

$$E(f)={\int }_{\mathcal{X}}c\left(x,f(x),\nabla f(x)\right)\mathrm{d}x,$$

$$v=v_1\wedge \dots \wedge v_k,$$

$$v=\sum _{\mathbf{i}\in I(d,k)}{v}^{\mathbf{i}}\cdot e_{{\mathbf{i}}_1}\wedge \dots \wedge e_{{\mathbf{i}}_k}=\sum _{\mathbf{i}\in I(d,k)}{v}^{\mathbf{i}}\cdot e_{\mathbf{i}},$$

$$\Vert v\Vert =\inf \left\{\sum _i|{\xi }_i|:{\xi }_i\text{ are simple},v=\sum _i{\xi }_i\right\}.$$

$$\Vert w{\Vert }^{\ast }=\sup \left\{⟨w,v⟩:v\text{ is simple },|v|\le 1\right\}.$$

$${\tau }_{{\mathcal{G}}_f}(z)=\frac{M(\nabla f({\pi }_{1}z))}{|M(\nabla f({\pi }_{1}z))|},$$

$$M(\xi )=(e_1+\xi e_1)\wedge \dots \wedge (e_n+\xi e_n),$$

$$v=\sum _{|\mathbf{i}|+|\mathbf{j}|=n}{v}^{\mathbf{i},\mathbf{j}}{e}_{\mathbf{i}}\wedge {\epsilon }_{\mathbf{j}},$$

$$\begin{array}{cccc}{v}^{\bar{0},0}& & & \\ \text{pagecolor}gray!10v^{1,1}& \text{pagecolor}gray!10v^{2,1}& & \\ \text{pagecolor}gray!10v^{1,2}& \text{pagecolor}gray!10v^{2,2}& v^{0,(1,2)}& \\ \text{pagecolor}gray!10v^{1,3}& \text{pagecolor}gray!10v^{2,3}& v^{0,(1,3)}& v^{0,(2,3)},\end{array}$$

$${\left[\xi (v)\right]}_{j,i}=(-1{)}^{n-i}{v}^{\bar{i},j}.$$

$${\Sigma }_1=\{v\in {\mathbf{\Lambda }}_n{\mathbf{R}}^{n+N}:v=M(\xi )\text{ for }\xi \in {\mathbf{R}}^{N\times n}\},$$

$$c(\xi )=\bar{c}(M(\xi ))\text{ for all }\xi \in {\mathbf{R}}^{N\times n}.$$

$$\Psi (z,v)=\{\begin{array}{ll}{v}^{\bar{0},0}\bar{c}({\pi }_{1}z,{\pi }_{2}z,v/v^{\bar{0},0}),& \text{if }{v}^{\bar{0},0}>0,\\ +\infty ,& \text{otherwise,}\end{array}$$

$${\int }_{\mathcal{X}}c(x,f(x),\nabla f(x))\mathrm{d}{\mathcal{L}}^n(x)$$

$$={\int }_{{\mathcal{G}}_f}\Psi (z,{\tau }_{{\mathcal{G}}_f}(z))\mathrm{d}{\mathcal{H}}^n(z),$$

$${\int }_{\mathcal{X}}c\left(x,f(x),\nabla f(x)\right)\mathrm{d}{\mathcal{L}}^n(x)$$

$$={\int }_{\mathcal{X}}\Psi \left(x,f(x),M(\nabla f(x))\right)\mathrm{d}{\mathcal{L}}^n(x)$$

$$={\int }_{{\mathcal{G}}_f}\Psi \left(z,M(\nabla f({\pi }_{1}z))\right)\frac{1}{|M(\nabla f({\pi }_{1}z))|}\mathrm{d}{\mathcal{H}}^n(z)$$

$$={\int }_{{\mathcal{G}}_f}\Psi \left(z,{\tau }_{{\mathcal{G}}_f}(z)\right)\mathrm{d}{\mathcal{H}}^n(z).$$

$${\int }_{\mathcal{M}}\omega :={\int }_{\mathcal{M}}⟨\omega (z),{\tau }_{\mathcal{M}}(z)⟩\mathrm{d}{\mathcal{H}}^k(z).$$

$$⟨d\omega (z),v_1\wedge \dots \wedge v_{k+1}⟩=\underset{h\to 0}{\lim }\frac{1}{h^{k+1}}{\int }_{\partial P}\omega ,$$

$${\int }_{\mathcal{M}}d\omega ={\int }_{\partial \mathcal{M}}\omega .$$

$$⟨{\pi }^{♯}\omega ,v_1\wedge ..\wedge v_k⟩=⟨\omega \circ \pi ,D_{v_1}\pi \wedge ..\wedge D_{v_k}\pi ⟩,$$

$$[[\mathcal{M}]](\omega )={\int }_{\mathcal{M}}\omega .$$

$$\partial T(\omega )=T(d\omega ),\text{ for all }\omega \in {\mathcal{D}}^{k-1}(U).$$

$$T(\omega )=0\text{ whenever }\mathrm{spt}(\omega )\subset V.$$

$${\pi }_{♯}T(\omega )=T({\pi }^{♯}\omega ),\text{ for all }\omega \in {\mathcal{D}}^k({\mathbf{R}}^q).$$

$$\mathbb{M}(T)=\sup \left\{T(\omega ):\omega \in {\mathcal{D}}^k(U),\Vert \omega (z){\Vert }^{\ast }\le 1\right\},$$

$$T(\omega )=\int ⟨\omega (z),\vec{T}(z)⟩\mathrm{d}\Vert T\Vert (z).$$

$$\mathbf{E}(T)=\int {\Psi }^{\ast \ast }\left({\pi }_{1}z,{\pi }_{2}z,\vec{T}(z)\right)\mathrm{d}\Vert T\Vert (z).$$

$$\mathbf{E}(T)=\underset{\omega \in \mathcal{K}}{\sup }T(\omega ),$$

$$\displaystyle\mathcal{K}=\Bigl{\{}$$

Full text

1 Introduction

We consider functionals of $C^{1}$-mappings $f:\mathcal{X}\to\mathcal{Y}$

[TABLE]

where $\mathcal{X}\subset\mathbf{R}^{n}$, $\mathcal{Y}\subset\mathbf{R}^{N}$ are bounded and open. The cost function $c\equiv c(x,y,\xi)$ is assumed to be a nonnegative (possibly nonconvex) continuous function on $\mathcal{X}\times\mathcal{Y}\times\mathbf{R}^{N\times n}$ that is polyconvex (see Def. 2) in the Jacobian matrix $\xi$.

This work is concerned with relaxation and global optimization of (1) when, both, dimension and codimension are possibly larger than one ($n>1$, $N>1$). This is expected to be difficult: In the discrete setting problems with $n=1$ or $N=1$ typically correspond to polynomial-time solvable shortest path ($n=1$) or graph cut ($N=1$) problems [cohen1997global, tsitsiklis, Ishikawa, schoenemann2010combinatorial], whereas for $n,N>1$, the arising multilabel problems with unordered label spaces are known to be NP-hard - see [li2016complexity]. Nevertheless, heuristic strategies have been shown to yield excellent results in tasks such as optical flow [fullflow] or shape matching [shekhovtsov2008efficient, chen2015robust]. In contrast to such well-established Markov random field (MRF) works [kolmogorov2006convergent, kolmogorov2007minimizing, kohli2008partial, shekhovtsov2008efficient, menze2015discrete, chen2015robust, fullflow, domokos2018mrf] we consider the way less explored continuous (infinite-dimensional) setting.

Our motivation partly stems from the fact that formulations in function space are very general and admit a variety of discretizations. Finite difference discretizations of continuous relaxations often lead to models that are reminiscent of MRFs [Zach-et-al-cvpr12], while piecewise-linear approximations are related to discrete-continuous MRFs [Zach-Kohli-eccv14], see [fix2014duality, moellenhoff-iccv-2017]. More recently, for the Kantorovich relaxation in optimal transport, approximations with deep neural networks were considered and achieved promising performance, for example in generative modeling [ACB17, seguy2017large].

We further argue that fractional (non-integer) solutions to a careful discretization of the continuous model can implicitly approximate an “integer” continuous solution. Therefore one can achieve accuracies that go substantially beyond the mesh size. The resulting models would be difficult to interpret and derive from a finite-dimensional viewpoint such that the continuous considerations are required for the final implementation. Also, formulations arising from continuous relaxations allow one to introduce isotropic smoothness potentials without reverting to higher-order terms in the cost, and, as we show in this work, one can impose general polyconvex regularizations using only local constraints. An example of a polyconvex function (which is in general nonconvex) is the surface area of the graph, sometimes referred to as “Beltrami regularization” in the image processing community, see e.g., [belt].

In contrast to the discrete multi-labeling setting, an important question is whether variational problems involving the energy (1) admit a minimizer. A fruitful approach to address this question is to suitably relax the notion of solution, thereby enlarging the search space of admissible candidates (“lifting the problem to a larger space”). The origins of this idea can be traced back111We refer the interested reader to the historical remarks in L. C. Young’s book on the calculus of variations [young2000lecture, pp. 122–123]. to the turn of the century, see Hilbert’s twentieth problem [hilbert]. An example of that principle is the celebrated Kantorovich relaxation [kantorovich1960] of Monge’s transportation problem [monge1781]. There, the search over maps $f:\mathcal{X}\to\mathcal{Y}$ is relaxed to one over probability measures on the product space $\mathcal{X}\times\mathcal{Y}$. Each map can be identified in that extended space with a measure concentrated on its graph. Existence of optimal transportation plans follows directly due to good compactness properties of the larger space. Furthermore, the nonlinearly constrained and nonconvex optimization problem is transformed into one of linear programming, leading to rich duality theories and fast numerical algorithms [PeCu18].

One may ask whether the relaxed solution in the extended space has certain regularity properties, for example whether it is the graph of a (sufficiently regular) map and thus can be considered a solution to the original (“unlifted”) problem. In the case of optimal transport, such regularity theory can be guaranteed under some assumptions [Vil08, San15]. Establishing existence and regularity for problems in which the cost additionally depends on the Jacobian (for example minimal surface problems) has been a driving factor in the development of geometric measure theory, see [morgan2016geometric] for an introduction. In this work, we will use ideas from geometric measure theory to pursue the above relaxation and lifting principle for the energy (1). The main idea is to reformulate the original variational problem as a shape optimization problem over oriented manifolds representing the graph of the map $f:\mathcal{X}\to\mathcal{Y}$ in the product space $\mathcal{X}\times\mathcal{Y}$ . To obtain a convex formulation we enlarge the search space from oriented manifolds to currents.

1.1 Related Work

A common strategy to solve problems involving (1) is to revert to local gradient descent minimization based on the Euler-Lagrange equations. But for nonconvex problems solutions might depend on the initialization and the computed stationary points may be quite suboptimal. Therefore, we pursue the aforementioned lifting of the energy (1) to currents. This lifting has been previously considered in geometric measure theory to establish the aforementioned existence and regularity theory for vectorial variational problems in a very broad setting, see e.g., [federer, federer1974real, aviles1991variational]. In contrast to such impressive theoretical achievements, this paper is concerned with a discretization and implementation.

There is also a variety of related applied works. The paper [windheuser2011geometrically] tackles the problem of bijective and smooth shape matching using linear programming. Similar to the present work, the authors also look for graph surfaces in $\mathcal{X}\times\mathcal{Y}$ but they consider the discrete setting and use a different notion of boundary operator. We study the continuous setting, but also our discrete formulation is quite different.

For $N=1$, the proposed continuous formulation specializes to [ABDM, PCBC-SIIMS]. To tackle the setting of $N>1$ in a memory efficient manner, Strekalovskiy et al. [Strekalovskiy-et-al-cvpr12, goldluecke2013tight, strekalovskiy-et-al-siims14] keep a collection of $N$ surfaces with codimension one under the factorization assumption that $\mathcal{Y}=\mathcal{Y}_{1}\times\ldots\times\mathcal{Y}_{N}$. In contrast, we consider only one surface of codimension $N$, we do not require an assumption on $\mathcal{Y}$, our approach is applicable to a larger class of functionals and we expect it to yield a tighter relaxation. The lifting approaches [lellmann-et-al-iccv2013, goldstein2012global] also tackle vectorial problems by considering the full product space, but are limited to total variation regularization (with the former allowing $\mathcal{Y}$ to be a manifold). The recent work [windheuser2016convex] is most related to the present one, however their relaxation considers a specific instance of (1). Moreover, the above works are based on finite difference discretizations of the continuous model. In contrast, the proposed discretization using discrete exterior calculus yields solutions beyond the mesh accuracy as in recent sublabel-accurate approaches. The latter are restricted to $N=1$ [moellenhoff-laude-cvpr-2016, moellenhoff-iccv-2017] or total variation regularization [laude16eccv]. Recent works also include extensions to total generalized variation or Laplacian regularization [strecke2018sublabel, vogt, loewe].

Recent approaches in shape analysis [solomon2016entropic, vestner2017product, vestner2017efficient] also operate in the product space $\mathcal{X}\times\mathcal{Y}$. However, these are based on local minimizations of the Gromov-Wasserstein distance [memoli2007use] and spectral variants thereof [memoli2009spectral] which leads to (nonconvex) quadratic assignment problems. While the goal to find a smooth (possibly bijective) map is similar, the formulations appear to be quite different. To alleviate the increased cost of the product space formulation, computationally efficient representations of densities in $\mathcal{X}\times\mathcal{Y}$ have been studied in the context of functional maps [ovsjanikov2012functional, rodola2018functional].

2 Notation and Preliminaries

Throughout this paper we will introduce notions from geometric measure theory, as they are not commonly used in the vision community. While the subject is rather technical, our aim is to keep the presentation light and to focus on the geometric intuition and aspects which are important for a practical implementation. We invite the reader to consult chapter 4 in the book [morgan2016geometric] and the chapter on exterior calculus in [crane2015discrete], which both contain many illuminating illustrations. For a more technical treatment we refer to [federer, KP08].

In the following, we denote a basis in $\mathbf{R}^{d}$ as $\{e_{1},\ldots,e_{d}\}$ with dual basis $\{\mathrm{d}x_{1},\ldots,\mathrm{d}x_{d}\}$ where $\mathrm{d}x_{i}:\mathbf{R}^{d}\to\mathbf{R}$ is the linear functional that maps every $x=(x_{1},\ldots,x_{d})$ to the $i$-th component $x_{i}$. Given an integer $k\leq d$, $I(d,k)$ are the multi-indices $\mathbf{i}=(i_{1},\ldots,i_{k})$ with $1\leq i_{1}<\ldots<i_{k}\leq d$.

As we will consider $n$-surfaces in $\mathcal{X}\times\mathcal{Y}\subset\mathbf{R}^{n+N}$, most of the time we set $d=n+N$ and $k=n$. To further simplify notation, we denote the basis vectors $\{e_{n+1},\ldots,e_{n+N}\}$ by $\{\varepsilon_{1},\ldots,\varepsilon_{N}\}$ and similarly refer to the dual basis as $\{\mathrm{d}x_{1},\ldots\mathrm{d}x_{n},\mathrm{d}y_{1},\ldots,\mathrm{d}y_{N}\}$. When it is clear from the context, we treat vectors $e_{i}\in\mathbf{R}^{n}$ and $\varepsilon_{i}\in\mathbf{R}^{N}$ in the sense that $e_{i}\simeq(e_{i},\mathbf{0}_{N})\in\mathbf{R}^{n+N}$, $\varepsilon_{i}\simeq(\mathbf{0}_{n},\varepsilon_{i})\in\mathbf{R}^{n+N}$. As an example, for $\nabla f(x)\in\mathbf{R}^{N\times n}$ we can define the expression $e_{i}+\nabla f(x)e_{i}$ and read it as $\left(e_{i},\nabla f(x)e_{i}\right)\in\mathbf{R}^{n+N}$.

2.1 Convex Analysis

The extended reals are denoted by $\overline{\mathbf{R}}=\mathbf{R}\cup\{+\infty\}$. For a finite-dimensional real vector space $V$ and $\Psi:V\to\overline{\mathbf{R}}$ we denote the convex conjugate as $\Psi^{*}:V^{*}\to\overline{\mathbf{R}}$ and the biconjugate as $\Psi^{**}:V\to\overline{\mathbf{R}}$. $\Psi^{**}$ is the largest lower-semicontinuous convex function below $\Psi$. In our notation, for functions with several arguments, the conjugate is always taken only in the last argument. As a general reference to convex analysis, we refer the reader to the books [hiriart2012fundamentals, Rockafellar:ConvexAnalysis].

2.2 Multilinear Algebra

The formalism of multi-vectors we introduce in this section is central to this work, as the idea of the relaxation is to represent the oriented graph of $f$ by a $k$-vectorfield (more precisely: a $k$-current) in the product space $\mathcal{X}\times\mathcal{Y}$. Basically, one can multiply $v_{i}\in\mathbf{R}^{d}$ to obtain an object

[TABLE]

called a simple $k$-vector in $\mathbf{R}^{d}$. The geometric intuition of simple $k$-vectors is, that they describe the $k$-dimensional space spanned by the $\{v_{i}\}$, together with an orientation and the area of the parallelotope given by the $\{v_{i}\}$. Thus, simple $k$-vectors can be thought of oriented parallelotopes as shown in orange in Fig. 1. In general, $k$-vectors are defined to be formal sums

[TABLE]

for coefficients $v^{\mathbf{i}}\in\mathbf{R}$. They form the vector space $\mathbf{\Lambda}_{k}\mathbf{R}^{d}$, which has dimension $\binom{d}{k}$.

The dual space $\mathbf{\Lambda}^{k}\mathbf{R}^{d}$ of $k$-covectors is defined analogously, with $\langle\mathrm{d}x_{\mathbf{i}},e_{\mathbf{j}}\rangle=\delta_{\mathbf{i}\mathbf{j}}$. We define for two $k$-vectors (and also for $k$-covectors) $v=\sum_{\mathbf{i}}v_{\mathbf{i}}e_{\mathbf{i}}$, $w=\sum_{\mathbf{i}}w_{\mathbf{i}}e_{\mathbf{i}}$ an inner product $\langle v,w\rangle=\sum_{\mathbf{i}}v_{\mathbf{i}}w_{\mathbf{i}}$ and norm $|v|=\sqrt{\langle v,v\rangle}$.

$k$-vectors (elements of $\mathbf{\Lambda}_{k}\mathbf{R}^{d}$) are called simple, if they can be written as the wedge product of $1$-vectors as in (2). Unfortunately, for $1<k<d-1$, not all $k$-vectors are simple and the set of simple $k$-vectors is a nonconvex cone in $\mathbf{\Lambda}_{k}\mathbf{R}^{d}$, called the Grassmann cone [busemann1963convex]. This is one aspect why the setting of $n>1$ and $N>1$ is more challenging.

Later on, we will consider a relaxation from the nonconvex set of simple $k$-vectors to general $k$-vectors. Naturally, for the relaxation to be good, we want the convex energy to be as large as possible on non-simple $k$-vectors. For the Euclidean norm, a good convex extension is the mass norm

[TABLE]

The dual norm is the comass norm given by:

[TABLE]

The mass norm can be understood as the largest norm that agrees with the Euclidean norm on simple $k$-vectors.

3 Lifting to Graphs in the Product Space

With the necessary preliminaries in mind, our goal is now to reparametrize the original energy (1) to the graph $\mathcal{G}_{f}\subset\mathcal{X}\times\mathcal{Y}$. As shown in Fig. 1, the graph is an oriented $n$-dimensional manifold in the product space with global parametrization $u(x)=(x,f(x))$.

Definition 1** (Orientation).**

If $\mathcal{M}\subset\mathbf{R}^{d}$ is a $k$-dimensional smooth manifold in $\mathbf{R}^{d}$ (possibly with boundary), an orientation of $\mathcal{M}$ is a continuous map $\tau_{\mathcal{M}}:\mathcal{M}\to\mathbf{\Lambda}_{k}\mathbf{R}^{d}$ such that $\tau_{\mathcal{M}}(z)$ is a simple $k$-vector with unit norm that spans the tangent space $T_{z}\mathcal{M}$ at every point $z\in\mathcal{M}$.

From differential geometry we know that the tangent space $T_{z}\mathcal{G}_{f}$ at $z=(x,f(x))$ is spanned by $\partial_{i}u(u^{-1}(z))=e_{i}+\nabla f(x)e_{i}$. Therefore, an orientation of $\mathcal{G}_{f}$ is given by

[TABLE]

where the map $M:\mathbf{R}^{N\times n}\to\mathbf{\Lambda}_{n}\mathbf{R}^{n+N}$ is given by

[TABLE]

and $\pi_{1}:\mathcal{X}\times\mathcal{Y}\to\mathcal{X}$ is the canonical projection onto the first argument. In order to derive the reparametrization, we have to connect a simple $n$-vector (representing an oriented tangent plane of the graph) with the Jacobian of the original energy. For that, we need an inverse of the map given in (7).

To derive such an inverse, we first introduce further helpful notations. For $\mathbf{i}\in I(m,l)$ we denote by $\bar{\mathbf{i}}\in I(m,m-l)$ the element which complements $\mathbf{i}$ in $\{1,2,\ldots,m\}$ in increasing order, denote $\bar{0}=\{1,\ldots,m\}$ and [math] as the empty multi-index. Every $v\in\mathbf{\Lambda}_{n}\mathbf{R}^{n+N}$ can be written as

[TABLE]

where $\mathbf{i}\in I(n,l)$, $\mathbf{j}\in I(N,l^{\prime})$, $l+l^{\prime}=n$. To give an example, the $\binom{5}{2}=10$ coefficients of a $2$-vector $v\in\mathbf{\Lambda}_{2}\mathbf{R}^{5}$ according to the notation (8) are:

[TABLE]

where we highlighted the $N\times n$ coefficients with $|\mathbf{j}|=1$. Now note that the vector $v=M(\xi)$ is by construction a simple $n$-vector with first component $v^{\bar{0},0}=1$. To any $v\in\mathbf{\Lambda}_{n}\mathbf{R}^{n+N}$ with $v^{\bar{0},0}=1$ we associate $\xi(v)\in\mathbf{R}^{N\times n}$ given by

[TABLE]

If and only if $v\in\mathbf{\Lambda}_{n}\mathbf{R}^{n+N}$ is simple with first component $v^{\bar{0},0}=1$ then $v=M(\xi(v))$. A proof is given in [GMS-CC, Vol. I, Ch. 2.1, Prop. 1]. Thus, on the set of simple $n$-vectors with first component $v^{\bar{0},0}=1$,

[TABLE]

the inverse of the map (7) is given by (10).

Using the above notations, we can define a generalized notion of convexity, which essentially states that there is a convex reformulation on $k$-vectors.

Definition 2** (Polyconvexity).**

A map $c:\mathbf{R}^{N\times n}\to\overline{\mathbf{R}}$ is polyconvex if there is a convex function $\bar{c}:\mathbf{\Lambda}_{n}\mathbf{R}^{n+N}\to\overline{\mathbf{R}}$ such that we have

[TABLE]

Equivalently one has that $c(\xi(v))=\bar{c}(v)$ for all $v\in\Sigma_{1}$. We also refer to the convex function $\bar{c}$ as a polyconvex extension.

In general, the polyconvex extension is not unique. Any convex function has an obvious polyconvex extension by (10), but as discussed in the previous section we would like the convex extension to be as large as possible for $v\notin\Sigma_{1}$. The largest polyconvex extension which agrees with the original function on $\Sigma_{1}$ can be formally defined using the convex biconjugate, but is often hard to explicitly compute. The mass norm (4) corresponds to such a construction.

Nevertheless, given any polyconvex extension, we can now reparametrize the original energy (1) on the oriented graph $\mathcal{G}_{f}$, as we show in the following central proposition.

Proposition 1**.**

Let $\bar{c}:\mathcal{X}\times\mathcal{Y}\times\mathbf{\Lambda}_{n}\mathbf{R}^{n+N}\to\overline{\mathbf{R}}$ be a polyconvex extension of the original cost $c$ in the last argument. Define the function $\Psi:\mathcal{X}\times\mathcal{Y}\times\mathbf{\Lambda}_{n}\mathbf{R}^{n+N}\to\overline{\mathbf{R}}$,

[TABLE]

where $\pi_{1}:\mathcal{X}\times\mathcal{Y}\to\mathcal{X}$ and $\pi_{2}:\mathcal{X}\times\mathcal{Y}\to\mathcal{Y}$ are the canonical projections onto the first and second argument. Then we can reparametrize (1) as follows:

[TABLE]

where the second integral is the standard Lebesgue integral with respect to the $n$-dimensional Hausdorff measure on $\mathbf{R}^{n+N}$ restricted to the graph $\mathcal{G}_{f}$.

Proof.

We directly calculate:

[TABLE]

The step from (15) to (16) uses that $\bar{c}$ is a polyconvex extension (so that we can apply (12)) and the fact that for $v=M(\nabla f(x))$ we have $v^{\bar{0},0}=1$. To arrive at (17), an application of the area formula [KP08, Corollary 5.1.13] suffices and for (18) we used positive one-homogenity of $\Psi$ and the definition of $\tau_{\mathcal{G}_{f}}$ in (6). ∎

Interestingly, the function (13) is convex and one-homogeneous in the last argument, as it is the perspective of a convex function. However, the search space of oriented graphs of $C^{1}$ mappings is nonconvex. Therefore we relax from oriented graphs to the larger set of currents, which we will introduce in the following section. Since currents form a vector space, we therefore obtain a convex functional over a convex domain.

4 From Oriented Graphs to Currents

Throughout this section, let $U\subset\mathbf{R}^{d}$ be an open set, which will later be a neighbourhood of $X\times Y\subset\mathbf{R}^{n+N}$, where $X=\operatorname{cl}(\mathcal{X})$, $Y=\operatorname{cl}(\mathcal{Y})$ are the closures of $\mathcal{X},\mathcal{Y}$. The main idea of our relaxation and the geometric intuitions of pushforward and boundary operator we introduce in this section are summarized in the following Fig. 2. Currents are defined in duality with differential forms, which we will briefly introduce in the following section.

4.1 Differential Forms

A differential form of order $k$ (short: $k$-form) is a map $\omega:U\to\mathbf{\Lambda}^{k}\mathbf{R}^{d}$. The support of a differential form $\operatorname{spt}\omega$ is defined as the closure of $\{z\in U:\omega(z)\neq 0\}$. Integration of a $k$-form over an oriented $k$-dimensional manifold is defined by

[TABLE]

A notion of derivative for $k$-forms is the exterior derivative $d\omega$, which is a $(k+1)$-form given by:

[TABLE]

where $\partial P$ is the oriented boundary of the parallelotope spanned by the $\{hv_{i}\}$ at point $z$.

To get an intuition, note that for $k=0$ this reduces to the familiar directional derivative $\langle d\omega(x),v_{1}\rangle=\lim_{h\to 0}\frac{1}{h}\left(\omega(x+hv_{1})-\omega(x)\right)$. With (19) and (20) in mind, one sees why Stokes’ theorem

[TABLE]

should hold intuitively. Given a map $\pi:\mathbf{R}^{d}\to\mathbf{R}^{q}$, the pullback $\pi^{\sharp}\omega$ of the $k$-form $\omega$ is determined by

[TABLE]

where $D_{v_{i}}\pi=\nabla\pi\cdot v_{i}$ and $\nabla\pi\in\mathbf{R}^{q\times d}$ is the Jacobian.

4.2 Currents

Denote the space of smooth $k$-forms with compact support on $U$ as $\mathcal{D}^{k}(U)$. Currents are elements of the dual space $\mathcal{D}_{k}(U)=\mathcal{D}^{k}(U)^{\prime}$, i.e., linear functionals acting on differential forms. As shown in Fig. 2(a), an oriented $k$-dimensional manifold $\mathcal{M}\subset U$ induces a current by

[TABLE]

However, since $\mathcal{D}_{k}(U)$ is a vector space, not all elements look like $k$-dimensional manifolds, see Fig. 2(d). The boundary of the $k$-current $T\in\mathcal{D}_{k}(U)$ is the $(k-1)$-current $\partial T\in\mathcal{D}_{k-1}(U)$ defined via the exterior derivative:

[TABLE]

Stokes’ theorem (21) ensures that for currents which are given by $k$-dimensional oriented manifolds, the boundary of the current agrees with the usual notion, see Fig. 2(b).

The support of a current, denoted by $\operatorname{spt}T$, is the complement of the biggest open set $V$ such that

[TABLE]

Given a map $\pi:\mathbf{R}^{d}\to\mathbf{R}^{q}$ the pushforward $\pi_{\sharp}T$ of the $k$-current $T\in\mathcal{D}_{k}(U)$ is given by

[TABLE]

Intuitively, it transforms the current using the map $\pi$, as illustrated in Fig. 2(a). The mass of a current $T\in\mathcal{D}_{k}(U)$ is

[TABLE]

and as expected $\mathbb{M}(\llbracket\mathcal{M}\rrbracket)=\mathcal{H}^{k}(\mathcal{M})$. We denote the space of $k$-currents with finite mass and compact support by $\mathbf{M}_{k}(U)$. These are representable by integration, meaning there is a measure $\|T\|$ on $U$ and a map $\vec{T}:U\to\mathbf{\Lambda}_{k}\mathbf{R}^{d}$ such that $\|\vec{T}(z)\|=1$ for $\|T\|$-almost all $z$ such that

[TABLE]

The decomposition (28) is crucial, and we will use it to define the relaxation in the next section.

4.3 The Relaxed Energy

We lift the original energy (1) to the space of finite mass currents $T\in\mathbf{M}_{n}(U)$ with $\operatorname{spt}T\subset X\times Y$ as follows:

[TABLE]

Since for $T=\llbracket\mathcal{G}_{f}\rrbracket$ we have $\vec{T}=\tau_{\mathcal{G}_{f}}$, \|T\|=\mathcal{H}^{n}\,\raisebox{-0.5468pt}{\reflectbox{\rotatebox[origin={br}]{-90.0}{\lnot}}}\,\mathcal{G}_{f} the desirable property $\mathbf{E}(\llbracket\mathcal{G}_{f}\rrbracket)=E(f)$ holds due to Prop. 1.

Note that in (29) we use the lower-semicontinuous regularization $\Psi^{**}$ which extends (13) at $v^{\bar{0},0}=0$ with the correct value. Interestingly, this point corresponds to the situation when the graph has vertical parts, which cannot occur for $C^{1}$ functions but can happen for general currents, see Fig. 2(c). In [Mora] it was shown that one can penalize such jumps in a way depending on the jump distance and direction. We will not consider such additional regularization due to space limitations, but remark that they could be integrated by adding further constraints to the following dual representation, which is a consequence of [GMS-CC, Vol. II, Sec. 1.3.1, Thm. 2].

Proposition 2**.**

For $T\in\mathbf{M}_{n}(U)$ with $\operatorname{spt}T\subset X\times Y$, we have the dual representation

[TABLE]

where the constraint is the closed and convex set

[TABLE]

The final relaxed optimization problem for (1) reads

[TABLE]

Depending on the kind of problem one wishes to solve, a different convex constraint set $\mathcal{C}$ should be considered. For example, in the case of variational problems with Dirichlet boundary conditions, we set

[TABLE]

where $S\in\mathbf{M}_{n-1}(U)$ is a given boundary datum. In case of Neumann boundary conditions, one constrains the support of the boundary to be zero inside the domain

[TABLE]

to exclude surfaces with holes, but allow the boundary to be freely chosen on $(\partial X)\times Y$. In case $n=N$, one can also consider the constraint set

[TABLE]

where the additional pushforward constraint encourages bijectivity. Notice also the similarity of (32) together with (35) to the Kantorovich relaxation in optimal transport.

Existence of minimizing currents to a similar problem as (32) in a certain space of currents (real flat chains) is shown in [federer1974real, §3.9]. For dimension $n=1$ or codimension $N=1$, the infimum is actually realized by a surface (integral flat chain) [federer1974real, §5.10, §5.12]. An adaptation of such theoretical considerations to our setting and conditions under which the relaxation is tight in the scenario $n>1$, $N>1$ is a major open challenge and left for future work.

5 Discrete Formulation

In this section we present an implementation of the continuous model (32) using discrete exterior calculus [Hirani2003]. We will base our discretization on cubes since they are easy to work with in high dimensions, but one could also use simplices. To define cubical meshes, we adopt some notations from computational homology [CH].

Definition 3** (Elementary interval and cube).**

An elementary interval is an interval $I\subset\mathbf{R}$ of the form $I=[l,l+1]$ or $I=\{l\}$ for $l\in\mathbf{Z}$. Intervals that consist of a single point are degenerate. An elementary cube is given by a product $\kappa=I_{1}\times\ldots\times I_{d}$, where each $I_{i}$ is an elementary interval. The set of elementary cubes in $\mathbf{R}^{d}$ is denoted by $K^{d}$.

For $\kappa\in K^{d}$, denote by $\dim\kappa\in\{1,\ldots,d\}$ the number of nondegenerate intervals. We denote $\mathbf{i}(\kappa)\in I(d,\dim\kappa)$ as the multi-index referencing the nondegenerate intervals.

Definition 4** (Cubical set).**

A set $Q\subset\mathbf{R}^{d}$ is a cubical set if it can be written as a finite union of elementary cubes.

Let $K_{k}^{d}(Q)=\{\kappa\in K^{d}:\kappa\subset Q,\dim\kappa=k\}$ be the set of $k$-dimensional cubes contained in $Q\subset\mathbf{R}^{d}$. A map $\phi:Q\to X\times Y$ will transform the cubical set to our domain. As we work with images, it will just be a mesh spacing, i.e., we set $\phi(z)=(h_{1}z_{1},\ldots,h_{d}z_{d})$.

Definition 5** ($k$-chains, $k$-cochains).**

We denote the space of finite formal sums of elements in $K_{k}^{d}(Q)$ with real coefficients as $\mathcal{C}_{k}(Q)$, called (real) $k$-chains. We denote the dual as $\mathcal{C}_{k}(Q)^{*}=\mathcal{C}^{k}(Q)$ and call the elements $k$-cochains.