# Lifting the Convex Conjugate in Lagrangian Relaxations:
A Tractable Approach for Continuous Markov Random Fields

Hartmut Bauermeister 44footnotemark: 4 equal contribution Emanuel
Laude11footnotemark: 1 Technical University of Munich, Department of
Informatics, Munich, Germany Thomas Möllenhoff RIKEN Center for Advanced
Intelligence Project, Tokyo, Japan.
Michael Moeller University of Siegen, Department of Electrical Engineering and
Computer Science, Siegen, Germany Daniel Cremers22footnotemark: 2

###### Abstract

Dual decomposition approaches in nonconvex optimization may suffer from a
duality gap. This poses a challenge when applying them directly to nonconvex
problems such as MAP-inference in a Markov random field (MRF) with continuous
state spaces. To eliminate such gaps, this paper considers a reformulation of
the original nonconvex task in the space of measures. This infinite-
dimensional reformulation is then approximated by a semi-infinite one, which
is obtained via a piecewise polynomial discretization in the dual. We provide
a geometric intuition behind the primal problem induced by the dual
discretization and draw connections to optimization over moment spaces. In
contrast to existing discretizations which suffer from a grid bias, we show
that a piecewise polynomial discretization better preserves the continuous
nature of our problem. Invoking results from optimal transport theory and
convex algebraic geometry we reduce the semi-infinite program to a finite one
and provide a practical implementation based on semidefinite programming. We
show, experimentally and in theory, that the approach successfully reduces the
duality gap. To showcase the scalability of our approach, we apply it to the
stereo matching problem between two images.

Keywords: Markov random fields $\cdot$ moment relaxation $\cdot$ sum of
squares $\cdot$ polynomial optimization $\cdot$ generalized conjugacy $\cdot$
optimal transport

AMS Subject Classification: 49N15 $\cdot$ 49M29 $\cdot$ 65K10 $\cdot$ 90C23
$\cdot$ 90C26 $\cdot$ 90C35

##  1 Introduction and Motivation

###  1.1 Lagrangian Relaxations

In this paper, our goal is to develop a convex optimization framework à la
dual decomposition for the MAP-inference problem in a continuous _Markov
random field_ (MRF). Continuous MRFs are versatile and therefore widely used
as a model in image processing, computer vision and machine learning [42, 53,
55, 10, 13, 50, 2, 51, 12, 3, 26, 36, 49]. Due to the continuous nature of the
state space, inference amounts to a continuous nonlinear optimization problem.
For an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ it takes the
following form:

$\displaystyle\min_{x\in\Omega^{\mathcal{V}}}\;F(x):=\sum_{u\in\mathcal{V}}f_{u}(x_{u})+\sum_{uv\in\mathcal{E}}f_{uv}(x_{u},x_{v}).$
(P)

The objective function $F$ is an additive composition of a separable part,
$\sum_{u\in\mathcal{V}}f_{u}(x_{u})$, with unary functions
${f_{u}:(\Omega\subset\mathbb{R}^{m})\to\mathbb{R}}$ and a coupling part
$\sum_{uv\in\mathcal{E}}f_{uv}(x_{u},x_{v})$ with pairwise functions
${f_{uv}:\Omega\times\Omega\to\mathbb{R}}$ and
$\emptyset\neq\Omega\subset\mathbb{R}^{m}$ compact.

The coupling term $\sum_{uv\in\mathcal{E}}f_{uv}(x_{u},x_{v})$ introduces a
major challenge for efficient optimization in particular when $\mathcal{V}$ is
large and therefore the problem is high-dimensional. For tractability, one
therefore seeks to find a decomposable reformulation of the problem. In
continuous optimization a viable approach is to derive a _Lagrangian
relaxation_ of the problem: The key idea is to introduce auxiliary variables
$x_{uv}\in\mathbb{R}^{m}\times\mathbb{R}^{m}$ for each edge $uv\in\mathcal{E}$
and linear constraints $x_{uv}=(x_{u},x_{v})$. Dualizing the linear
constraints with Lagrange multipliers
$\lambda_{uv}\in\mathbb{R}^{m}\times\mathbb{R}^{m}$ one arrives at the
Lagrangian dual problem which falls within the regime of convex optimization
algorithms that can exploit its separable dual structure. Examples include
subgradient ascent or the _primal-dual hybrid gradient_ (PDHG) method [9].

However, the approach suffers from potentially large duality gaps since
problem (P) is nonconvex in general. Indeed, it can be shown that a direct
Lagrangian relaxation of (P) leads to the following “naive” convexification of
the original problem:

$\displaystyle\min_{x\in(\mathbb{R}^{m})^{\mathcal{V}}}\sum_{u\in\mathcal{V}}f_{u}^{**}(x_{u})+\sum_{uv\in\mathcal{E}}f_{uv}^{**}(x_{u},x_{v}).$
(1)

With some abuse of notation, here, $f_{u}$ and $f_{uv}$ attain the value
$+\infty$ whenever $x_{u}\not\in\Omega$ or $(x_{u},x_{v})\not\in\Omega^{2}$.
Then $f_{u}^{**}$ are the convex biconjugates which correspond to the largest
convex lower semicontinuous under-approximations to $f_{u}$. Such component-
wise convex envelopes can produce inaccurate or even trivial convex under-
approximations to the global objective $F$.

###  1.2 The Local Marginal Polytope Relaxation

To overcome such duality gaps, we propose to reformulate the original
nonlinear problem in terms of an infinite-dimensional linear program over the
space of Radon measures and apply the Lagrangian relaxation afterwards.

The following key lemma reveals, that every minimization problem can be
equivalently formulated in terms of an infinite-dimensional linear program. In
the following, we denote $\mathcal{P}(X)$ as the space of Borel probability
measures on $\emptyset\neq X\subset\mathbb{R}^{m}$ compact, and we denote by
$\langle\mu,f\rangle=\int_{X}f(x)\,\mathrm{d}\mu(x)$.

######  Lemma 1.1.

Let $f:X\to\mathbb{R}$ be lower semicontinuous (lsc) with $\emptyset\neq
X\subset\mathbb{R}^{m}$ compact. Then we have

$\displaystyle\min_{x\in
X}f(x)=\min_{\mu\in\mathcal{P}(X)}~{}\langle\mu,f\rangle,$ (2)

and $x^{*}\in\operatorname*{arg\,min}_{x\in X}f(x)$ is a solution to
$\min_{x\in X}f(x)$ if and only if $\delta_{x^{*}}$ is a minimizer of
$\min_{\mu\in\mathcal{P}(X)}\langle\mu,f\rangle$.

###### Proof.

The result follows immediately via compactness of $X$ and the properties of
probability measures: By compactness of $X$ and since $f$ is lsc relative to
$X$ we know $x^{*}$ exists. Let $f_{\min}=f(x^{*})$. By the properties of the
Lebesgue integral we have:

$\int_{X}f(x)\,\mathrm{d}\mu(x)\geq\int_{X}f_{\min}\,\mathrm{d}\mu(x)=f_{\min}\int_{X}1\,\mathrm{d}\mu(x)=\min_{x\in
X}f(x),$

for all $\mu\in\mathcal{P}(X)$ and
$\min_{x\in\mathbb{R}^{n}}f(x)=\langle\delta_{x^{*}},f\rangle=f(x^{*})$. ∎

Let $\pi_{u}:\Omega^{\mathcal{V}}\to\Omega$,
$\pi_{uv}:\Omega^{\mathcal{E}}\to\Omega^{2}$ denote the canonical projections
onto the $u\textsuperscript{th}$ resp. $u\textsuperscript{th}$ and
$v\textsuperscript{th}$ components and $\pi_{u}\sharp\mu$ is the pushforward
of $\mu$ w.r.t. $\pi_{u}$ defined by:
$(\pi_{u}\sharp\mu)(A)=\mu((\pi_{u})^{-1}(A))$ for all $A\subset\Omega$. Then,
in our case, applying the reformulation from Lem. 1.1 directly to the cost
function $F$ with $X=\Omega^{\mathcal{V}}$ we obtain by linearity:

$\displaystyle\min_{\mu\in\mathcal{P}(\Omega^{\mathcal{V}})}\;\langle\mu,F\rangle$
$\displaystyle=\min_{\mu\in\mathcal{P}(\Omega^{\mathcal{V}})}\;\left\langle\mu,\sum_{u\in\mathcal{V}}f_{u}\circ\pi_{u}+\sum_{uv\in\mathcal{E}}f_{uv}\circ\pi_{uv}\right\rangle$
$\displaystyle=\min_{\mu\in\mathcal{P}(\Omega^{\mathcal{V}})}\;\sum_{u\in\mathcal{V}}\langle\pi_{u}\sharp\mu,f_{u}\rangle+\sum_{uv\in\mathcal{E}}\langle\pi_{uv}\sharp\mu,f_{uv}\rangle.$

This relaxation is known as the _full marginal polytope relaxation_ which is,
however, intractable if $\mathcal{V}$ is large as one minimizes over
probability measures on the product space $\mathcal{P}(\Omega^{\mathcal{V}})$.
Instead, we consider the following linear programming relaxation of (P) which
is also referred to as the _local marginal polytope relaxation_ [26, 14, 48,
35] which is more tractable as the optimization variable lies in the product
space of probability measures $\mathcal{P}(\Omega)^{\mathcal{V}}$:

$\displaystyle\inf_{\mu\in\mathcal{P}(\Omega)^{\mathcal{V}}}\mathcal{F}(\mu):=\sum_{u\in\mathcal{V}}\langle\mu_{u},f_{u}\rangle+\sum_{uv\in\mathcal{E}}\operatorname{OT}_{f_{uv}}(\mu_{u},\mu_{v}).$
(R-P)

Here, $\operatorname{OT}_{f_{uv}}$ denotes the optimal transportation [15]
with marginals $(\mu_{u},\mu_{v})$ and cost $f_{uv}$ defined by

$\displaystyle\operatorname{OT}_{f_{uv}}(\mu_{u},\mu_{v})=\inf_{\mu_{uv}\in\Pi(\mu_{u},\mu_{v})}\langle
f_{uv},\mu_{uv}\rangle.$ (3)

The constraint set $\Pi(\mu_{u},\mu_{v})$ consists of all Borel probability
measures on $\Omega^{2}$ with specified marginals $\mu_{u}$ and $\mu_{v}$:

$\displaystyle\Pi(\mu_{u},\mu_{v})=\left\\{\mu_{uv}\in\mathcal{P}(\Omega^{2}):\pi_{u}\sharp\mu_{uv}=\mu_{u},\;\pi_{v}\sharp\mu_{uv}=\mu_{v}\right\\}.$
(4)

Note that for finite state-spaces (i.e., $|\Omega|<\infty$), the set
$\mathcal{P}(\Omega)$ can be identified with the standard
$(|\Omega|-1)$-dimensional probability simplex and $\Pi(\mu_{u},\mu_{v})$ with
the set of nonnegative $|\Omega|\times|\Omega|$ matrices whose rows and
columns sum up to $\mu_{v}$ and $\mu_{u}$. In that case, the linear program
given in (R-P) is equivalent to the well-known finite-dimensional local
marginal polytope relaxation for MRFs, which is for example studied in [52].
We further remark that the case of finite $\Omega$ has been extensively
studied in the literature, see, e.g., [16, 47] for recent overviews.

For the more challenging setting of continuous state-spaces, a major
difficulty stems from the fact that the linear programming relaxation (R-P) is
an infinite-dimensional optimization problem posed in the space of Borel
probability measures. Perhaps due to this difficulty, discrete MRF approaches
are still routinely applied despite the continuous nature of $\Omega$. This is
typically done by considering a finite sample approximation of $\Omega$, so
that the infinite-dimensional linear program reduces to a finite-dimensional
one. This, however, may lead to discretization errors and comes with an
exponential complexity in the dimension of $\Omega$.

Due to the fact that
$\Pi(\delta_{x},\delta_{x^{\prime}})=\\{\delta_{(x,x^{\prime})}\\}$ one sees
that restricting $\mu_{u}$ (and therefore $\mu_{uv}$) to be Dirac probability
measures, the formulation (R-P) reduces to the original problem (P). As one
instead considers the larger convex set of all probability measures, it is a
_relaxation_ which lower bounds (P), i.e., we have the following important
relation:

$\eqref{eq:opt_relax}\leq\eqref{eq:mrf}.$ (5)

When considering a relaxation, the immediate question arises whether this
lower bound is attained, i.e., the relaxation is tight and therefore the above
inequality holds with equality. For finite and ordered $\Omega$ the situation
is well-understood, see [52]. For infinite $\Omega\subset\mathbb{R}^{m}$, to
our knowledge, the situation is less clear, unless $f_{u}$, $f_{uv}$ are
convex. An exception is, when $f_{u}$, $f_{uv}$ are continuous piecewise
linear: Then the continuous formulation is equivalent to a discrete MRF, see
[14, Thm. 6], and existing results from the discrete case carry over. Despite
the lack of theoretical tightness guarantees, in experiments, the infinite-
dimensional local marginal polytope relaxation produces accurate lower
relaxations, in particular when $\Omega$ is univariate.

A Lagrangian relaxation to the infinite-dimensional problem is obtained by
dualizing the marginalization constraints (4) for each $e\in\mathcal{E}$ with
Lagrange multipliers $\lambda_{e}\in\mathcal{C}(\Omega)^{2}$. This is
equivalent to a substitution of the optimal transportation with its dual
formulation. Adopting the approach of [14] a finite-dimensional problem is
obtained by approximating the Lagrange multipliers in terms of finite linear
combinations of certain basis functions. These approximations are chosen in
such a way that the classical Lagrangian relaxation and the discrete approach,
described above, are special cases of the considered framework: Let $D$ be the
Lagrange dual function. Then for any hierarchy of increasingly expressive dual
subspaces
$\Lambda_{1}\subset\Lambda_{2}\subset\cdots\subset\mathcal{C}(\Omega)$, where
$\Lambda_{1}$ is the space of affine functions, the induced hierarchy of dual
problems

$\displaystyle\eqref{eq:classical_lagr_relaxation}$
$\displaystyle=\max_{\lambda\in((\Lambda_{1})^{2})^{\mathcal{E}}}\;D(\lambda)\leq\max_{\lambda\in((\Lambda_{2})^{2})^{\mathcal{E}}}\;D(\lambda)\leq\cdots\leq\max_{\lambda\in(\mathcal{C}(\Omega)^{2})^{\mathcal{E}}}\;D(\lambda)=\eqref{eq:opt_relax}\leq\eqref{eq:mrf},$

leads to a reduction of the duality gap. For a piecewise polynomial hierarchy
with increasing degrees and/or number of pieces we will show that the gap
eventually vanishes as the number of pieces and/or the degree goes to
$\infty$.

###  1.3 Contributions and Overview

Expanding upon existing related work our contributions are threefold:

_Theory_ : We derive and study the primal optimization problem corresponding
to the dual discretization and draw connections to optimization over moment
spaces and aspects from variational analysis. In the same way a convex duality
follows from a linear separation argument, a nonconvex duality can be
formulated in terms of nonlinear separation boundaries. After a “lifting” to a
higher-dimensional space these nonlinear separation boundaries become linear.
This is analogous to the idea of feature maps for linear classifiers and gives
rise to our notion of a “lifted” convex conjugate. More formally, we show that
under a certain extremality condition, satisfied by a polynomial
discretization with degree at least $2$, the generalized biconjugate of a
potentially discontinuous, nonconvex function equals (up to closure) the
original function. As a consequence, under extremality, the discretized energy
agrees with the original energy, when restricted to a certain continuum of
extreme points associated to Diracs. And, for every point in the continuous
domain of the original energy such a correspondence exists. In contrast, in a
piecewise linear discretization there only exist finitely many extreme points
which makes such a one-to-one correspondence impossible. In addition, for a
spline-based dual discretization we show that the nonconvex primal-dual gap
vanishes as the degree and/or number of pieces goes to infinity.

_First tractable formulation and implementation_ : We exploit optimal
transport duality theory, and the aforementioned duality relations to
optimization over moment spaces to derive, to our knowledge, the first
tractable polynomial discretization of the infinite-dimensional local marginal
polytope relaxation. More precisely, after discretization the problem amounts
to a semi-infinite program which can be transformed into a separable finite-
dimensional semidefinite program applying concepts from algebraic geometry,
such as nonnegativity certificates of polynomials. This allows us to derive an
efficient first-order primal-dual algorithm which, due to the separable
problem structure, can be parallelized on a GPU to handle large problems.

_Experimental evaluation and large scale capability_ : Experimentally, we
provide numerical evidence for the strict reduction of the duality-gap and
demonstrate the advantages of the nonlinear approximation over existing
liftings in the literature in Sec. 5. We implement our algorithm on a GPU and
apply it to the nonconvex problem of large-scale stereo matching. The results
show that increasing the dual subspace improves both, the dual energy, and the
quality of a recovered primal solution from a discretized measure.

#### Notation

For a compact nonempty set $\Omega\subset\mathbb{R}^{m}$, denote by
$\mathcal{P}(\Omega)$ the space of Borel probability measures on $\Omega$ and
by $\mathcal{M}(\Omega)$ the space of Radon measures on $\Omega$. The convex
cone of nonnegative Radon measures is denoted by $\mathcal{M}_{+}(\Omega)$.
Furthermore let $\mathcal{C}(\Omega)$ be the space of continuous functions on
$\Omega$. We will write $\langle
f,\mu\rangle=\int_{\Omega}f(x)\,\mathrm{d}\mu(x)$, for
$f\in\mathcal{C}(\Omega)$ and $\mu\in\mathcal{M}(\Omega)$. For a set
$C\subset\mathbb{R}^{m}$ we denote by
$\iota_{C}:\mathbb{R}^{m}\to\overline{\mathbb{R}}=\mathbb{R}\cup\\{\pm\infty\\}$
the indicator function with $\iota_{C}(x)=0$ if $x\in C$ and
$\iota_{C}(x)=\infty$ if $x\notin C$ and by $\sigma_{C}(x)=\sup_{y\in
C}\langle x,y\rangle$ the support function of $C$ at $x\in\mathbb{R}^{m}$. We
denote by $C^{*}=\\{y\in\mathbb{R}^{m}:\langle y,x\rangle\geq 0,\forall x\in
C\\}$ the dual cone of $C$. These notions are defined analogously for the
topologically paired spaces $\mathcal{M}(\Omega)$ and $\mathcal{C}(\Omega)$.
The convex hull $\operatorname{con}C$ of a set $C\subset\mathbb{R}^{m}$ is the
smallest convex set that contains $C$. Equivalently, $\operatorname{con}C$ is
the set of all finite convex combinations of points in $C$. With some abuse of
notation we write $\operatorname{con}f$ for a function $f$ to denote the
largest convex function below $f$. Let $\delta_{x}$ denote the Dirac measure
centered at $x\in\mathbb{R}^{m}$. We write lsc for lower semicontinuous and
denote the epigraphical closure of some function $f$ by $\operatorname{cl}f$.
Let $\llbracket\cdot\rrbracket$ denote the Iverson bracket, where $\llbracket
P\rrbracket=1$ if $P$ is true and $\llbracket P\rrbracket=0$, otherwise. For
some extended real-valued function $f:\mathbb{R}^{m}\to\overline{\mathbb{R}}$
let $f^{*}(y)=\sup_{x\in\mathbb{R}^{m}}\langle x,y\rangle-f(x)$ denote the
Fenchel conjugate of $f$ at $y$ and $f^{**}=(f^{*})^{*}$ the Fenchel
biconjugate.

###  1.4 Related Work

#### MRFs

Fix and Agarwal [14] are the first to propose a dual subspace approximation of
the infinite-dimensional local marginal polytope relaxation. This allows for a
generalization and unified treatment of convex piecewise linear formulations
for MRFs with continuous state spaces proposed in [55, 54]. However, for the
polynomial case, due to the semi-infinite problem structure, the existence of
an implementable algorithm beyond subgradient ascent is posed as an open
question by [14]. As a consequence, no numerical results are presented. Dual
subspace approximations for related models in a spatially continuous setting
have been considered by [23] and further employed for the global optimization
of vector- or manifold-valued optimization problems in imaging and vision [24,
45, 46]. In these papers, piecewise linear dual approximations are used and
the techniques developed in the present work allow one to go beyond the
piecewise linear case in a tractable way. Besides applications in imaging and
vision, MRFs with continuous state spaces have for example been used in
protein folding [26].

#### Lifting

Lifting to measures for global constrained polynomial optimization is used in
[18, 19] to reformulate the infinite-dimensional _linear program_ (LP) in
terms of a linear objective over the _semidefinite programming_ (SDP)
characterization of the finite-dimensional space of moments. Its SDP-dual is
connected to nonnegativity certificates of polynomials [17, 41, 39]. As a key
difference to [18, 19] we consider possibly nonpolynomial objective functions
which results in a generally nonlinear formulation over the space of moments.
Similar techniques have been considered recently in [20, 11] to estimate
Lipschitz constants in neural networks.

#### Lagrangian relaxation and decoupling by lower relaxations

A component-wise lifting in problems with a partially separable structure
results in a generalized dual decomposition approach. Dual decomposition and
Lagrangian relaxation are general principles in optimization and appear as a
useful tool across many disciplines, see [4, Ch. 5–6], [22] and the references
therein for an overview. Traditionally, decomposition methods are applied
directly in the nonconvex setting without lifting, unlike the generalized
scheme we present in this paper which is based on a lifted reformulation.
Without lifting this typically results in a component-wise convex lower
envelope. More recently, [40] consider a homotopy method based on component-
wise Lasry–Lions envelopes which specialises to a component-wise convex
envelope in a certain limit case.

#### Generalized conjugacy and duality

Our notion of a lifted convex conjugate is closely related to generalized
conjugate functions originally due to [25]. It was utilized by [1] to study
nonconvex dualities, expanding upon the work of [32] on augmented Lagrangians
for nonconvex optimization. The specific case of quadratic conjugate functions
was developed by [27] as an analytic tool in a seminal proof that the proximal
subgradient map of a lower semicontinuous extended real-valued function is
monotone if and only if the function is convex. Our motivation, shared with
[32], is to eliminate duality gaps in nonconvex optimization. We offer a
somewhat complementary view through a primal-type lifting approach which
allows us to establish the connection to optimization in spaces of measures.

##  2 Dual Discretization for the Continuous MRF

###  2.1 A Reduced Dual Formulation

In contrast to previous approaches [14], we restrict ourselves to metric
pairwise terms $f_{uv}(x,y)=d(x,y)$ which eventually leads to a different dual
formulation and turns out more tractable: Then, the optimal transportation
$\operatorname{OT}_{f_{uv}}(\mu_{u},\mu_{v})$ in Problem (R-P) is the
Wasserstein-$1$ distance $W_{1}^{d}(\mu_{u},\mu_{v})$ induced by the metric
$d$ between $\mu_{u}$ and $\mu_{v}$. Furthermore we assume that
$f_{u}:\Omega\to\mathbb{R}$ is lsc and $\Omega\subset\mathbb{R}^{m}$ is a
compact nonempty set. In particular this implies that $f_{u}$ is proper and
bounded from below.

Thanks to optimal transport duality theory [44] we are therefore able to
obtain a more compact dual formulation, which is instrumental to derive a
tractable implementation for a piecewise polynomial discretization later on:
More precisely, we substitute
$\operatorname{OT}_{f_{uv}}(\mu_{u},\mu_{v})=W_{1}^{d}(\mu_{u},\mu_{v})$ in
Problem (R-P) with its dual formulation

$\displaystyle
W_{1}^{d}(\mu_{u},\mu_{v})=\sup_{\lambda\in\mathrm{Lip}_{d}(\Omega)}\int\lambda(x)\,\mathrm{d}(\mu_{u}-\mu_{v})(x),$
(6)

where $\lambda$ is $1$-Lipschitz with respect to the metric $d$, i.e.,

$\displaystyle\lambda\in\mathrm{Lip}_{d}(\Omega)=\\{\lambda:\Omega\to\mathbb{R}:|\lambda(x)-\lambda(y)|\leq
d(x,y)\\}.$ (7)

In contrast to the general dual formulation of optimal transport, which
involves two Lagrange multipliers per edge, interacting via the constraint
set, the Wasserstein-$1$ dual involves only a single dual variable for each
edge, that satisfies a Lipschitz constraint. We show in Sec. 4 that even
though this formulation still involves infinitely many constraints, thanks to
convex algebraic geometry [5], there exists a tractable finite representation
of this constraint set in terms of SDP.

For notational convenience, we assign an arbitrary orientation to the edges in
$\mathcal{G}$. After introducing a graph divergence operator
$\operatorname{Div}:\mathcal{C}(\Omega)^{\mathcal{E}}\to\mathcal{C}(\Omega)^{\mathcal{V}}$
defined by:

$-(\operatorname{Div}\lambda)_{u}=\sum_{v:(u,v)\in\mathcal{E}}\lambda_{(u,v)}-\sum_{v:(v,u)\in\mathcal{E}}\lambda_{(v,u)},$
(8)

an interchange of $\min$ and $\sup$ yields the following reduced dual problem,
which is the starting point for further discussion and a tractable
implementation:

$\displaystyle\sup_{\lambda\in\mathcal{C}(\Omega)^{\mathcal{E}}}D(\lambda):=-\sum_{u\in\mathcal{V}}\sigma_{\mathcal{P}(\Omega)}(-f_{u}+(\operatorname{Div}\lambda)_{u})-\sum_{e\in\mathcal{E}}\iota_{\mathcal{K}}(\lambda_{e}),$
(R-D)

where $\mathcal{K}=\mathrm{Lip}_{d}(\Omega)$ and
$\sigma_{\mathcal{P}(\Omega)}(-f_{u}+(\operatorname{Div}\lambda)_{u})=\sup_{\mu\in\mathcal{P}(\Omega)}\langle\mu,-f_{u}+(\operatorname{Div}\lambda)_{u}\rangle$
is the support function of $\mathcal{P}(\Omega)$ at
$-f_{u}+(\operatorname{Div}\lambda)_{u}$, which, thanks to Lem. 1.1, can be
rewritten:

$\displaystyle-\sigma_{\mathcal{P}(\Omega)}(-f_{u}+(\operatorname{Div}\lambda)_{u})=\min_{x\in\Omega}(f_{u}-(\operatorname{Div}\lambda)_{u})(x).$
(9)

Under mild assumptions we have strong duality:

######  Proposition 2.1.

Let $\emptyset\neq\Omega\subset\mathbb{R}^{m}$ be compact, let
$f_{u}:\Omega\to\mathbb{R}$ be lsc, and $f_{uv}=d$, where
$d:\Omega^{2}\to\mathbb{R}$ is lsc and a metric. Then, the following strong
duality holds:

$\displaystyle\eqref{eq:opt_relax}=\eqref{eq:dualMPTsmall},$ (10)

and a maximizer of (R-D) exists.

###### Proof.

Define $f:\mathcal{M}(\Omega)^{\mathcal{V}}\to\overline{\mathbb{R}}$

$f(\mu):=\sum_{u\in\mathcal{V}}\langle
f_{u},\mu_{u}\rangle+\iota_{\mathcal{P}(\Omega)}(\mu_{u}),$

which is convex, proper and lsc due to [37, Lem. 1.1.3]. Likewise, define the
functional $g:\mathcal{M}(\Omega)^{\mathcal{E}}\to\overline{\mathbb{R}}$

$g(\nu):=\sum_{e\in\mathcal{E}}\sup_{\lambda_{e}\in\mathcal{K}}\langle\lambda_{e},\nu_{e}\rangle=\sum_{e\in\mathcal{E}}\sigma_{\mathcal{K}}(\nu_{e}),$

which is proper convex lsc, as it is a pointwise supremum over linear
functionals. Define
$\nabla:\mathcal{M}(\Omega)^{\mathcal{V}}\to\mathcal{M}(\Omega)^{\mathcal{E}}$

$(\nabla\mu)_{(u,v)}=\mu_{u}-\mu_{v},$

which is bounded and linear. Then we compute the convex conjugates
$f^{*}:\mathcal{C}(\Omega)^{\mathcal{V}}\to\overline{\mathbb{R}}$ to

$f^{*}(\theta)=\sum_{u\in\mathcal{V}}\sigma_{\mathcal{P}(\Omega)}(\theta_{u}-f_{u}),$

and the functional
$g^{*}:\mathcal{C}(\Omega)^{\mathcal{E}}\to\overline{\mathbb{R}}$

$g^{*}(\lambda)=\sum_{e\in\mathcal{E}}\iota_{\mathcal{K}}(\lambda_{e}),$

and
$\nabla^{*}:\mathcal{C}(\Omega)^{\mathcal{E}}\to\mathcal{C}(\Omega)^{\mathcal{V}}=-\operatorname{Div}$.

Choose $x\in\Omega$ and define $\mu:=(\delta_{x})^{\mathcal{V}}$. Then
$f(\mu)=\sum_{u\in\mathcal{V}}f_{u}(x)<\infty$. In addition we have
$\langle\lambda_{e},\delta_{x}-\delta_{x}\rangle=0$ for all
$\lambda_{e}\in\mathcal{C}(\Omega)$ and therefore $g(\nabla\mu)=0$. Now
consider a weakly∗-convergent sequence
$\mathcal{M}(\Omega)^{\mathcal{E}}\ni\nu^{t}\overset{\ast}{\rightharpoonup}\nabla\mu$.
This means for all $e\in\mathcal{E}$ and $\lambda_{e}\in\mathcal{C}(\Omega)$
we have
$\langle\lambda_{e},\nu_{e}^{t}\rangle\to\langle\lambda_{e},(\nabla\mu)_{e}\rangle=0$.
In particular this implies that $g(\nu^{t})\to 0$, and therefore $g$ is
continuous at $\nabla\mu$. Then we can invoke the Fenchel–Rockafellar duality
Theorem [30, Thm. 1] and obtain that
$\eqref{eq:opt_relax}=\eqref{eq:dualMPTsmall}$ and a maximizer of (R-D)
exists. ∎

As discussed in Sec. 1 the next step in our strategy to obtain a tractable
formulation is to restrict $\lambda_{(u,v)}\in\mathcal{C}(\Omega)$ in Problem
(R-D) to a subspace $\Lambda=\langle\varphi_{0},\ldots,\varphi_{n}\rangle$
spanned by basis functions $\varphi_{k}\in\mathcal{C}(\Omega)$ and instead
consider $\sup_{\lambda\in\Lambda^{\mathcal{E}}}D(\lambda)$.

Before we derive the primal problem corresponding to the discretized dual
problem we discuss the duality between finite-dimensional subspaces of
$\mathcal{C}(\Omega)$ and certain equivalence classes of measures and in
particular moment spaces in the next section.

###  2.2 Duality of $\Lambda\subset\mathcal{C}(\Omega)$

Due to Riesz’ theorem the dual space of continuous functions on a compact set
$\mathcal{C}(\Omega)$ is given by the space of Radon measures
$\mathcal{M}(\Omega)$. Note that conversely for any $\lambda\in\Lambda$ we
obtain a linear functional on $\mathcal{M}(\Omega)$ by mapping
$\mu\in\mathcal{M}(\Omega)$ to
$\mu\mapsto\langle\mu,\lambda\rangle=\int\lambda(x)\,\mathrm{d}\mu(x)$.
Although any $\mu\in\mathcal{M}(\Omega)$ maps to a linear functional on
$\Lambda$, this mapping is not injective. In particular there exist distinct
measures $\mu\neq\nu$ such that
$\langle\mu,\lambda\rangle=\langle\nu,\lambda\rangle$ for every
$\lambda\in\Lambda$, i.e., $\mu$ and $\nu$ induce the same linear functional
on $\Lambda$. Instead, invoking [34, Thm. 4.9] the dual space $\Lambda^{*}$
can be related to a quotient space of $\mathcal{M}(\Omega)$ as follows:
$\Lambda^{*}$ is isometrically isomorphic to a quotient space on
$\mathcal{M}(\Omega)$ via

$\displaystyle\Lambda^{*}\cong\mathcal{M}(\Omega)/\Lambda^{0}:=\left\\{\mu+\Lambda^{0}\colon\mu\in\mathcal{M}(\Omega)\right\\},$
(11)

where
$\Lambda^{0}:=\left\\{\mu\in\mathcal{M}\colon\langle\mu,\lambda\rangle=0\text{
for all $\lambda\in\Lambda$}\right\\}$ is the annihilator of $\Lambda$. Any
$\mu\in\mathcal{M}(\Omega)$ generates an equivalence class in
$\mathcal{M}/\Lambda^{0}$ denoted by $[\mu]$ and the corresponding element in
$\Lambda^{*}$ is denoted by $[\mu]_{\Lambda}$.

Further restricting to the cone of functions in $\Lambda$ that are nonnegative
on $\Omega$, the corresponding dual cone has a nice interpretation: It can be
identified with a certain space of moments: Given a measure
$\mu\in\mathcal{M}_{+}(\Omega)$, for a particular choice of basis functions
$\\{\varphi_{0},\varphi_{1},\ldots,\varphi_{n}\\}$, we refer to
$\langle\mu,\varphi_{k}\rangle=\int\varphi_{k}(x)\,\mathrm{d}\mu(x)$ as the
$k$th _moment_ of $\mu$ and likewise the moment space is

$\displaystyle(\mathcal{M}_{\Lambda})_{+}=\left\\{y\in\mathbb{R}^{n+1}:\exists\mu\in\mathcal{M}_{+}(\Omega):y_{k}=\langle\mu,\varphi_{k}\rangle\right\\}\subset\mathbb{R}^{n+1}.$
(12)

This is discussed in Sec. 4 in the context of implementation. The above result
then shows that the moment space $(\mathcal{M}_{\Lambda})_{+}$ identifies the
space of equivalence classes of nonnegative measures, where two nonnegative
measures are considered equivalent if they have the same moments.

For the remainder of this paper, for simplicity, we fix the basis
$\\{\varphi_{0},\varphi_{1},\ldots,\varphi_{n}\\}$ and identify the dual
variable $\lambda_{(u,v)}$ with its vector of coefficients
$p_{(u,v)}\in\mathbb{R}^{n+1}$.

###  2.3 Dual Discretization for the Metric Pairwise MRF

In view of the Wasserstein-$1$ dual in Eq. (6) constant components in
$\lambda_{(u,v)}$ cancel out and therefore, assuming $\varphi_{0}$ is the
constant function, we can choose $p_{0}=0$. As a consequence we can represent
$(\operatorname{Div}\lambda)_{u}$ in terms of its coefficients
$(\operatorname{Div}_{\Lambda}p)_{u}$ w.r.t. the basis functions
$\\{\varphi_{1},\ldots,\varphi_{n}\\}$, for an appropriately chosen linear
mapping
$\operatorname{Div}_{\Lambda}:(\mathbb{R}^{n})^{\mathcal{E}}\to(\mathbb{R}^{n})^{\mathcal{V}}$.

We identify the basis functions $\varphi_{k}$ as the component functions of
the mapping $\varphi:\Omega\to\mathbb{R}^{n}$. Then, substituting the
representation $\lambda_{(u,v)}=\langle p_{(u,v)},\varphi(\cdot)\rangle$ in
the support function we obtain thanks to Eq. (9) the following discretized
dual problem:

$\sup_{p\in(\mathbb{R}^{n})^{\mathcal{E}}}D_{\Lambda}(p):=\sum_{u\in\mathcal{V}}~{}\min_{x\in\Omega}\;-\big{\langle}(\operatorname{Div}_{\Lambda}p)_{u},\varphi(x)\big{\rangle}+f_{u}(x)-\sum_{e\in\mathcal{E}}\iota_{\mathcal{K}_{\Lambda}}(p_{e}),$
(dR-D)

where the constraint set $\mathcal{K}_{\Lambda}$ is the set of coefficients of
Lipschitz functions in $\Lambda$:

$\displaystyle\mathcal{K}_{\Lambda}:=\\{p\in\mathbb{R}^{n}:\langle
p,\varphi(\cdot)\rangle\in\mathrm{Lip}_{d}(\Omega)\\}.$ (13)

In order to derive the primal problem note that for any
$f_{u}:\Omega\to\mathbb{R}$ the expression
$\sup_{x\in\Omega}\langle(\operatorname{Div}_{\Lambda}p)_{u},\varphi(x)\rangle-
f_{u}(x)$, (up to the presence of $\varphi(x)$) resembles the form of a convex
conjugate. Indeed, exploiting the notion of an extended real-valued function
this is the convex conjugate of a lifted version
$f_{\Lambda}:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ of
$f:\Omega\to\mathbb{R}$ which is defined by

$\displaystyle f_{\Lambda}(y)=\begin{cases}f(x),&\text{if $y=\varphi(x)$ for
some $x\in\Omega$,}\\\ \infty,&\text{otherwise},\end{cases}$ (14)

and gives rise to our notion of a _lifted convex conjugate_ which is studied
in depth in Sec. 3.

Since the “lifting” to a higher dimensional space happens through the
application of the mapping $\varphi$, whose component functions are the
elements of the basis, it is intuitive to refer to $\varphi$ as the lifting
aka feature map in this context. Likewise we refer to the curve described by
$x\mapsto\varphi(x)$ as the moment curve.

With this construction at hand and via an application of Fenchel–Rockafellar
duality we are ready to state the primal problem of the discretized dual
problem in terms of the lifted biconjugates $(f_{u})_{\Lambda}^{**}$:

$\displaystyle\min_{y\in(\mathbb{R}^{n})^{\mathcal{V}}}~{}\mathcal{F}_{\Lambda}(y):=\sum_{u\in\mathcal{V}}(f_{u})_{\Lambda}^{**}(y_{u})+\sum_{e\in\mathcal{E}}\sigma_{\mathcal{K}_{\Lambda}}\big{(}(\nabla_{\Lambda}y)_{e}\big{)},$
(dR-P)

where $\nabla_{\Lambda}:=-\operatorname{Div}_{\Lambda}^{*}$ is the negative
adjoint of the graph divergence, the graph gradient operator. A comparison of
(dR-P) with the convex relaxation (1) shows the effect of the dual
discretization, where the classical convex biconjugates are replaced with
biconjugates of the “lifted” functions $(f_{u})_{\Lambda}$. Indeed, without
lifting, i.e., $\varphi(x)=x$, we recover the classical biconjugates and
therefore the convex relaxation (1).

In experiments, we will discretize the dual problem with a piecewise
polynomial family of functions. For intervals, the following proposition shows
that either by increasing the number of pieces or the degree of the polynomial
the primal-dual gap can be reduced. In our case we approximate the Lipschitz
dual variable in terms of a Lipschitz spline. As a consequence existing
results such as [14, Thm. 2] do not apply. Instead, we use a construction
based on Bernstein-polynomials. Then the result follows from [6, Thm. 1].

######  Proposition 2.2.

Assume that $\Omega=[a,b]\subset\mathbb{R}$, $a<b$, and let the metric $d$ be
given by $d(x,y)=|x-y|$. Furthermore, let $\Lambda\subset\mathcal{C}(\Omega)$
be the space spanned by continuous piecewise polynomials on intervals
$[t_{i},t_{i+1}]$ defined by a regularly spaced grid with nodes given by
$t_{i}=a+(b-a)\cdot(i-1)/K$, $i=1,\ldots,K+1$. Then the optimality gap
satisfies:

$\eqref{eq:opt_relax}-\eqref{eq:discrete_dual}=\mathcal{O}(1/(K\cdot\sqrt{\deg})),$

where $\deg$ is the degree of the polynomial on each piece.

###### Proof.

We consider the discretized dual problem (dR-D) where $K$ is the set of
coefficients corresponding to $1$-Lipschitz piecewise polynomials on $[a,b]$
of degree $\deg$ with $k$ pieces. Also recall that we have the following
relations between the dual and primal problems:
$\eqref{eq:discrete_primal}=\eqref{eq:discrete_dual}\leq\eqref{eq:opt_relax}=\eqref{eq:dualMPTsmall}$.

Now, let us denote a maximizer of (R-D) as
$\lambda^{*}\in\mathrm{Lip}_{d}(\Omega)^{\mathcal{E}}$. Existence of such a
dual maximizer follows by Prop. 2.1. Then, one has for any
$\lambda\in\Lambda^{\mathcal{E}}$:

$\displaystyle\min_{x\in\Omega}~{}f_{u}(x)-(\operatorname{Div}\lambda^{*})_{u}(x)$
$\displaystyle=\min_{x\in\Omega}~{}f_{u}(x)-(\operatorname{Div}\lambda)_{u}(x)-(\operatorname{Div}\lambda^{*})_{u}(x)+(\operatorname{Div}\lambda)_{u}(x)$
$\displaystyle\leq\min_{x\in\Omega}~{}f_{u}(x)-(\operatorname{Div}\lambda)_{u}(x)+\|-(\operatorname{Div}(\lambda^{*}-\lambda))_{u}\|_{\infty}.$
(15)

This allows us to bound the optimality gap by:

$\displaystyle\eqref{eq:dualMPTsmall}-\eqref{eq:discrete_dual}$
$\displaystyle\leq\sum_{u\in\mathcal{V}}\|-(\operatorname{Div}(\lambda-\lambda^{*}))_{u}\|_{\infty}$
$\displaystyle\leq\sum_{u\in\mathcal{V}}|d(u)|\cdot\sup_{e\in\mathcal{E}}~{}\|\lambda_{e}-\lambda_{e}^{*}\|_{\infty}$
$\displaystyle\leq
2|\mathcal{E}|\cdot\sup_{e\in\mathcal{E}}\|\lambda_{e}-\lambda_{e}^{*}\|_{\infty},$
(16)

where $d(u)$ denotes the degree of the vertex $u$.

For a $L$-Lipschitz function $f:[0,1]\to\mathbb{R}$ there exists a Bernstein
polynomial $p:[0,1]\to\mathbb{R}$ with $p(0)=f(0)$ and $p(1)=f(1)$ such that
$\|p-f\|_{\infty}\leq\frac{3L}{2}\deg^{-1/2}$ [8, Thm. 2.6]. By [6, Thm. 1],
this polynomial is $L$-Lipschitz as well. For each $e\in\mathcal{E}$ we pick
the coefficients of the function $\lambda_{e}$ such that it approximates the
optimal dual variable $\lambda_{e}^{*}$ with such a polynomial individually on
each interval $[t_{i},t_{i+1}]$. Then one obtains an overall $1$-Lipschitz
polynomial with the following bound:

$\|\lambda_{e}-\lambda_{e}^{*}\|_{\infty}\leq\frac{3(b-a)}{2k\sqrt{\deg}}.$
(17)

Inserting this into (16) yields via
$\eqref{eq:opt_relax}=\eqref{eq:dualMPTsmall}$:

$\eqref{eq:opt_relax}-\eqref{eq:discrete_dual}\leq
3|\mathcal{E}|\frac{(b-a)}{k\sqrt{\deg}},$ (18)

which gives the stated $\mathcal{O}(1/(k\cdot\sqrt{\deg}))$ rate. ∎

##  3 Lifted Convex Conjugates

###  3.1 Lifting and Moments

(a) Standard Lifting

(b) Sublabel Lifting [21]

(c) Polynomial Lifting (Ours)

Figure 1: Different finite-dimensional approximations $\mathcal{P}_{\Lambda}$
of the infinite-dimensional space of probability measures
$\mathcal{P}([-1,1])$ with $\Omega=[-1,1]$. Left and middle: $2$-dimensional
probability simplex and right: Monomial moment space
$\operatorname{con}\\{(x,x^{2},x^{3}):x\in[-1,1]\\}$ of degree $3$. The
approximations are obtained as the convex hulls of the black curves
$x\mapsto\varphi(x)$ for 3 different choices of $\varphi$. From left to right
cf. Ex. 3.3, Ex. 3.4 and Ex. 3.5. In all cases, the black curves themselves
correspond to Dirac measures and the convex hulls of the curves correspond to
the space of probability measures. In contrast to the simplex that only has a
finite number of extreme points, the monomial moment curve comprises a
continuum of extreme points so that no Dirac measure on the monomial moment
curve can be expressed as a convex combination of other Diracs.

A fundamental question in the study of the lifted biconjugate
$(f_{u})_{\Lambda}^{**}$ is the charcterization of its domain
$\operatorname{dom}(f_{u})_{\Lambda}^{**}=\operatorname{con}\varphi(\Omega)$,
which, assuming $f_{u}$ is bounded from below, equals the convex hull of the
graph of $\varphi$. For that purpose we define the probability moment space,
as the set of all vectors $y\in\mathbb{R}^{n}$ for which there exists a
probability measure $\mu\in\mathcal{P}(\Omega)$ such that the $k^{\text{th}}$
component of $y$ is the $k^{\text{th}}$ moment of $\mu$:

$\displaystyle\mathcal{P}_{\Lambda}=\left\\{y\in\mathbb{R}^{n}:\exists\mu\in\mathcal{P}(\Omega),y_{k}=\langle\mu,\varphi_{k}\rangle\right\\},$
(19)

In other words $\mathcal{P}_{\Lambda}$ is the set of all “infinite” convex
combinations of points $y\in\varphi(\Omega)$ that belong to the graph of
$\varphi$. The following relation to $(\mathcal{M}_{\Lambda})_{+}$ holds: If
$\varphi_{0}\equiv 1$, a vector
$(y_{0},y_{1},\ldots,y_{n})\in(\mathcal{M}_{\Lambda})_{+}$ if and only if
$(y_{1},\ldots,y_{n})\in\mathcal{P}_{\Lambda}$. Also see Lem. 4.1.

Furthermore, note that for some $x\in\Omega$ the moment vector
$\varphi(x)\\!=\int\varphi(x^{\prime})\,\mathrm{d}\delta_{x}(x^{\prime})$ of a
Dirac measure $\delta_{x}$ sets up a correspondence between lifted points
$\varphi(x)$ and the Dirac measure $\delta_{x}$ iteself. Thus the moment space
can be written equivalently as the set of all finite convex combinations of
moment vectors of Dirac measures, in the same way the unit simplex is the
convex hull of the unit vectors:

######  Proposition 3.1.

Let $\emptyset\neq\Omega$ be compact and
$\Lambda=\langle\varphi_{0},\ldots,\varphi_{n}\rangle\subset\mathcal{C}(\Omega)$
with $\varphi_{0}\equiv 1$. Then every moment vector of a probability measure
is a finite convex combination of moment vectors of Dirac measures, i.e.

$\displaystyle\operatorname{con}({\varphi(\Omega)})=\mathcal{P}_{\Lambda},$
(20)

and $\mathcal{P}_{\Lambda}$ is nonempty and compact.

###### Proof.

Note that $\mathcal{P}_{\Lambda}$ is convex and bounded. It is also closed: To
this end consider a sequence $y^{t}\to y$ with
$y^{t}\in\mathcal{P}_{\Lambda}$. This means for any $t$ there exists
$\mu^{t}\in\mathcal{P}(\Omega)$ with
$y_{k}^{t}=\int_{\Omega}\varphi_{k}(x)\,\mathrm{d}\mu^{t}(x)$ for
$k=1,\ldots,n$. Since $\Omega$ is compact, due to Prokhorov’s theorem, see
[37, Sec. 1.1], there exists a weakly∗ convergent subsequence
$\mu^{t_{j}}\overset{\ast}{\rightharpoonup}\mu$ and, hence,
$\int_{\Omega}\varphi_{k}(x)\,\mathrm{d}\mu^{t_{j}}(x)=y_{k}^{t_{j}}\to\int_{\Omega}\varphi_{k}(x)\,\mathrm{d}\mu(x)=y_{k}$.
Therefore $y\in\mathcal{P}_{\Lambda}$. Next we show identity of the support
functions of the convex sets $\operatorname{con}({\varphi(\Omega)})$ and
$\mathcal{P}_{\Lambda}$ as this implies the equality of the sets. To this end
let $f\in\Lambda=\langle\varphi_{0},\ldots,\varphi_{n}\rangle$ with
$\varphi_{0}=1$. We write $f(x)=\langle\varphi(x),f\rangle$. Assume that
$f_{0}=0$. We have the identities:

$\displaystyle\sigma_{\varphi(\Omega)}(-f)=-\min_{x\in\Omega}~{}\langle\varphi(x),f\rangle=-\min_{\mu\in\mathcal{P}(\Omega)}~{}\int_{\Omega}f(x)\,\mathrm{d}\mu(x)=-\min_{y\in\mathcal{P}_{\Lambda}}~{}\langle
y,f\rangle=\sigma_{\mathcal{P}_{\Lambda}}(-f),$

where the first equality follows by the definition of the support function,
the second equality by Lem. 1.1 and the third equality by the definition of
$\mathcal{P}_{\Lambda}$, the identity

$\min_{\mu\in\mathcal{P}(\Omega)}~{}\int_{\Omega}f(x)\,\mathrm{d}\mu(x)=\min_{\mu\in\mathcal{P}(\Omega)}\sum_{k=1}^{n}f_{k}\int_{\Omega}\varphi_{k}(x)\,\mathrm{d}\mu(x),$

and the substitution $y_{k}=\int_{\Omega}\varphi_{k}(x)\,\mathrm{d}\mu(x)$.
Since for each such $f$ the support functions are equal we have equality of
the support functions of $\varphi(\Omega)$ and $\mathcal{P}_{\Lambda}$. Since
$\Omega$ is compact and $\varphi$ continuous and the convex hull of a compact
set stays compact, cf. [33, Cor. 2.30], we can replace $\varphi(\Omega)$ with
its convex hull $\operatorname{con}\varphi(\Omega)$ and the conclusion
follows. ∎

The proof of the above proposition also reveals that for $f\in\Lambda$, the
lifted biconjugate is a linear function over the moment space,

$\displaystyle f_{\Lambda}^{**}(y)=\langle
y,f\rangle+\iota_{\mathcal{P}_{\Lambda}}(y),$ (21)

whose minimization is actually equivalent to minimizing the original function.

Specializing $\varphi_{k}$ to the monomial basis and $\Omega$ to a set defined
via polynomial inequalities, this yields the formulation proposed by [18, 19]
for constrained polynomial optimization.

More generally, a meaningful notion of moments $y_{k}$ is induced by a lifting
map $\varphi$ which is a homeomorphism between $\Omega$ and $\varphi(\Omega)$.

######  Definition 3.2 (lifting map).

Let $\Omega$ be compact and be nonempty. Then we say the mapping
$\varphi:\Omega\to\mathbb{R}^{n}$ is a lifting map if $\varphi$ is continuous
on $\Omega$ and injective with continuous inverse
$\varphi^{-1}:\varphi(\Omega)\to\mathbb{R}^{m}$.

It is instructive to discuss possible choices for $\varphi$ including the ones
that correspond to existing discretizations for the continuous MRF such as the
discrete approach and the piecewise linear approach. In the latter two cases,
the moment space is merely the unit simplex. In that sense, the monomial
moment space can be interpreted as a “nonlinear probability simplex”, which,
in contrast to the unit simplex, has infinitely many extreme points, see Fig.
1. For the same reason, as we will see in the course of this section, it
better suits the continuous nature of our optimization problem.

The discrete sampling-based approach is recovered by the following choice of
$\varphi$:

######  Example 3.3.

We discretize the interval $\Omega=[a,b]$, $a<b$ and re-define
$\Omega:=\\{t_{1},\ldots,t_{n}\\}$ with $t_{k}\in[a,b]$, $t_{k}<t_{k+1}$. For
any $t_{k}\in\Omega$ let $\varphi(t_{k})=e_{k}$ with $e_{k}\in\mathbb{R}^{n}$
being the $k^{\text{th}}$ unit vector. As a result $\varphi$ is the canonical
basis and spans the space of discrete functions
$f:\\{t_{1},\ldots,t_{n}\\}\to\mathbb{R}$ and the moment space
$\operatorname{con}({\varphi(\Omega)})=\mathcal{P}_{\Lambda}$ is given by the
unit simplex.

The above example can be extended by assigning any points
$x\in(t_{k},t_{k+1})$ to points on the connecting line between two
corresponding Diracs which results in a more continuous formulation:

######  Example 3.4.

Let $\Omega=[a,b]$, $a<b$: Let $t_{k}<t_{k+1}$ and $t_{1}=a$, $t_{n}=b$ be a
sequence of knots that subdivide the interval $\Omega$ into $n-1$ subintervals
$[t_{k},t_{k+1}]=:\Omega_{k}$. We define $\varphi(x)=\alpha
e_{k}+(1-\alpha)e_{k+1}$, with $e_{k}\in\mathbb{R}^{n}$ being the
$k^{\text{th}}$ unit vector, $\alpha\in[0,1]$ such that $\alpha
t_{k}+(1-\alpha)t_{k+1}=x$. This yields a so-called $2$-sparse lifting map
that has been used in the related work of [21]. It spans the space $\Lambda$
of univariate piecewise linear functions on $\Omega$ and the moment space
$\operatorname{con}({\varphi(\Omega)})=\mathcal{P}_{\Lambda}$ is given by the
unit simplex.

For $\varphi_{k}$ being chosen as the monomials we obtain the classical notion
of moments:

######  Example 3.5.

For $\varphi_{0}=1$, the space of univariate polynomials
$\Lambda=\mathbb{R}[x]$ with maximum degree $n$ is spanned by the monomials:

$\displaystyle\varphi(x)=(x,x^{2},\ldots,x^{n}),$ (22)

and $\operatorname{con}({\varphi(\Omega)})=\mathcal{P}_{\Lambda}$ is the
monomial moment space.

######  Example 3.6.

Identifying $\Omega=\\{x\in\mathbb{R}^{2}\colon\lVert x\rVert_{2}=1\\}$ with
the complex unit circle $\Omega\cong\\{z\in\mathbb{C}\colon\lvert
z\rvert=1\\}$ and again assuming $\varphi_{0}=1$, the mapping

$\displaystyle\varphi(z)=(\operatorname{Re}(z),\operatorname{Im}(z),\cdots,\operatorname{Re}(z^{n}),\operatorname{Im}(z^{n}))\in\mathbb{R}^{2n}$
(23)

spans the space $\Lambda$ of real trigonometric polynomials of maximum degree
$n$. Parametrizing elements $z\in\Omega$ via the bijection between
$z=e^{i\omega}$ and its angle $\omega\in[0,2\pi)$ the components of $\
\varphi$ are the Fourier basis functions, which define the Carathéodory curve
[7]. The convex hull
$\operatorname{con}({\varphi(\Omega)})=\mathcal{P}_{\Lambda}$ is the
trigonometric moment space.

###  3.2 Extremal Subspaces

In contrast to the piecewise linear lifting, the polynomial lifting is
extremal in the sense that no Dirac measure on the monomial moment curve can
be expressed as a convex combination of other Diracs, which sets up an one-to-
one correspondece between Diracs $\delta_{x}$ and lifted points $\varphi(x)$.
This is illustrated in Fig. 1.

More formally, we call a lifting map $\varphi$ an extremal curve if each
$y\in\varphi(\Omega)$ is an extreme point of
$\operatorname{con}(\varphi(\Omega))$ defined according to [31, Sec. 18].

######  Definition 3.7 (extreme points).

Let $C$ be a convex set and $y\in C$. Then $y$ is called an extreme point of
$C$ if there is no way to express $y$ as a convex combination
$y=(1-\alpha)x+\alpha z$ of $x,z\in C$ and $0<\alpha<1$, except by taking
$y=x=z$.

######  Definition 3.8 (extremal moment curve).

Let $\Omega$ be compact and be nonempty. Then we say the mapping
$\varphi:\Omega\to\mathbb{R}^{n}$ is an extremal moment curve if $\varphi$ is
a lifting map and any point $y\in\varphi(\Omega)\subset\mathbb{R}^{n}$ is an
extreme point of $\operatorname{con}\varphi(\Omega)$.

Via a change of basis it becomes clear that the definition of extremality is
independent of a specific choice of a basis for $\Lambda$. Therefore
extremality is rather a property of the subspace $\Lambda$. This also
motivates the following lemma which shows that extremality is inherited along
a hierarchy $\Theta\subset\Lambda\subset\mathcal{C}(\Omega)$.

######  Lemma 3.9 (extremal subspaces).

Let $\emptyset\neq\Omega\subset\mathbb{R}^{m}$ be compact. Let
$\Theta\subset\Lambda\subset\mathcal{C}(\Omega)$ be a hierarchy of finite-
dimensional subspaces of the space of continuous functions
$\mathcal{C}(\Omega)$. Let $\Theta=\langle\theta_{1},\ldots,\theta_{n}\rangle$
such that $\theta:\Omega\to\mathbb{R}^{n}$ is an extremal curve. Then
$\Lambda$ is spanned by an extremal curve as well.

###### Proof.

$\\{\theta_{1},\ldots,\theta_{n}\\}$ is a basis of $\Theta\subset\Lambda$ and
therefore linearly independent. Since $\Lambda$ is finite-dimensional in view
of the basis extension theorem $\\{\theta_{1},\ldots,\theta_{n}\\}$ can be
extended to a basis
$\varphi=(\theta_{1},\ldots,\theta_{n},\psi_{1},\ldots,\psi_{k})$ of $\Lambda$
with vectors $\psi_{i}\in G$, where $\operatorname{span}G=\Lambda$ and
$|G|<\infty$ such that $\operatorname{span}\varphi=\Lambda$.

Now choose $y\in\Omega$ and consider $\varphi(y)$. Let $\alpha\in(0,1)$ and
$\varphi(y)=\alpha x+(1-\alpha)z$ for
$x,z\in\operatorname{con}\varphi(\Omega)$. Then there exist coefficients
$\alpha_{i}$, $\beta_{k}>0$ such that
$x=\sum_{i=1}^{m}\alpha_{i}\varphi(x^{(i)})$ and
$z=\sum_{k=1}^{l}\beta_{k}\varphi(z^{(k)})$, $z^{(k)},x^{(i)}\in\Omega$ with
$\sum_{i=1}^{m}\alpha_{i}=1$, $\sum_{k=1}^{l}\beta_{k}=1$. This implies that
$\theta(y)=\alpha\sum_{i=1}^{m}\alpha_{i}\theta(x^{(i)})+(1-\alpha)\sum_{k=1}^{l}\beta_{k}\theta(z^{(k)})$.
Extremality of $\theta$ implies that
$\theta(y)=\sum_{i=1}^{l}\alpha_{i}\theta(x^{(i)})=\sum_{k=1}^{l}\beta_{k}\theta(z^{(k)})$
and therefore $\theta(y)=\theta(x^{(i)})=\theta(z^{(k)})$. Since $\theta$ is
an extremal curve it is injective and therefore $y=z^{(k)}=x^{(i)}$. This
implies that $\varphi(y)=x=z$. ∎

######  Lemma 3.10 (extremality of quadratic subspace).

Let $\emptyset\neq\Omega\subset\mathbb{R}^{m}$ be compact. Let
$\varphi(x)=(x_{1},x_{2},\ldots,x_{m},\|x\|^{2})$. Assume that
$\\{x\mapsto\langle
a,\varphi(x)\rangle:a\in\mathbb{R}^{m+1}\\}\subset\Lambda\subset\mathcal{C}(\Omega)$.
Then $\Lambda$ is spanned by an extremal curve.

###### Proof.

Now choose $y\in\Omega$. Let $\alpha\in(0,1)$ and $\varphi(y)=\alpha
x+(1-\alpha)z$ for $x,z\in\operatorname{con}\varphi(\Omega)$. Then there exist
coefficients $\alpha_{i}$, $\beta_{k}>0$ such that
$x=\sum_{i=1}^{m}\alpha_{i}\varphi(x^{(i)})$ and
$z=\sum_{k=1}^{l}\beta_{k}\varphi(z^{(k)})$, $z^{(k)},x^{(i)}\in\Omega$ with
$\sum_{i=1}^{m}\alpha_{i}=1$, $\sum_{k=1}^{l}\beta_{k}=1$.

Now choose $f(x)=\|x-y\|^{2}=\|x\|^{2}-2\langle y,x\rangle+\|y\|^{2}$. Then we
have for $a=(-2y,1)$ and $a_{0}=\|y\|^{2}$:

$\displaystyle 0=f(y)=\langle\varphi(y),a\rangle+a_{0}$
$\displaystyle=\left\langle\alpha\sum_{i=1}^{m}\alpha_{i}\varphi(x^{(i)})+(1-\alpha)\sum_{k=1}^{l}\beta_{k}\varphi(z^{(k)}),a\right\rangle+a_{0}$
$\displaystyle=\sum_{i=1}^{m}\alpha\cdot\alpha_{i}\cdot\langle\varphi(x^{(i)}),a\rangle+\sum_{k=1}^{l}(1-\alpha)\cdot\beta_{k}\cdot\langle\varphi(z^{(k)}),a\rangle+a_{0}$
$\displaystyle=\sum_{i=1}^{m}\alpha\cdot\alpha_{i}\cdot
f(x^{(i)})+\sum_{k=1}^{l}(1-\alpha)\cdot\beta_{k}\cdot f(z^{(k)})$

As $f(x^{(i)})>0$ for $x^{(i)}\neq y$ and $\alpha>0$ it holds that
$x^{(i)}\neq y$ implies $\alpha_{i}=0$. The same is true for $z^{(k)}$ and
$\beta_{k}$. Hence $x=\varphi(y)=z$, and therefore
$\varphi:\Omega\to\mathbb{R}^{m+1}$ is an extremal curve. In view of Lem. 3.9
$\Lambda$ is spanned by an extremal curve. ∎

This shows that in a piecewise polynomial discretization with degree at least
$2$ the corresponding basis inherits the extreme point property from the
extremality of the subspace of quadratic functions.

(a) Standard Lifting

(b) Sublabel Lifting [21]

(c) Polynomial Lifting (Ours)

Figure 2: Lifted version $f_{\Lambda}$ of the function $f$ for 3 different
choices of $\varphi$. In the top row the same nonconvex function $f$ (blue
curves) is depicted. The colored surfaces in the bottom row correspond to the
lifted biconjugates $f_{\Lambda}^{**}$, where the gray shadow areas correspond
to their domains $\operatorname{dom}f_{\Lambda}^{**}=\mathcal{P}_{\Lambda}$
and the blue curves correspond to the lifted cost $f_{\Lambda}$, see Eq. (14).
The black curves in the bottom row correspond to the moment curve described by
Diracs $\varphi(\Omega)$, i.e., the domain of $f_{\Lambda}$. From left to
right, cf. Ex. 3.3, Ex. 3.4 and Ex. 3.5. The nonlinear supporting dual
functions $q_{\lambda,\beta}(x)=\langle\varphi(x),\lambda\rangle-\beta$ (red
curves) to $f$ in the top row (middle and right), are transformed into linear
supporting hyperplanes $l_{\lambda,\beta}(y)=\langle y,\lambda\rangle-\beta$
(red surfaces) to $f_{\Lambda}^{**}$ in the bottom row through the feature map
$\varphi$. In the language of kernel methods, these functions can be
interpreted as the nonlinear decision boundaries that separate individual
points on the graph of $f$ from the epigraph of $f$. Only in the most right
case such a separation is possible. As a result the polynomial lifting
preserves the nonconvex cost function $f_{\Lambda}^{**}\circ\varphi=f$ on
$\varphi(\Omega)$ whereas the $2$-sparse lifting (middle) only leads to a
piecewise convex under-approximation $f_{\Lambda}^{**}\circ\varphi\leq f$.

Extremal curves are key to preserve the cost function when restricted to the
set of discretized Diracs $\varphi(x)$:

######  Theorem 3.11.

Let $\Omega\subset\mathbb{R}^{m}$ be nonempty and compact and let
$f:\Omega\to\mathbb{R}$ be lsc. Furthermore, let
$\varphi:\Omega\to\mathbb{R}^{n}$ be an extremal curve. Then we have

$\displaystyle f_{\Lambda}^{**}\circ\varphi=f$ (24)

on $\Omega$. In addition, we have that
$\operatorname{con}f_{\Lambda}=f_{\Lambda}^{**}$.

###### Proof.

Since $f$ is lsc and $\Omega$ compact $f$ is bounded from below, i.e., there
is $\gamma>-\infty$ so that $f(x)\geq\gamma$ for all $x\in\Omega$. We have
$\operatorname{dom}f_{\Lambda}=\varphi(\Omega)\subset\mathbb{R}^{n}$ and in
view of [33, Prop. 2.31] it holds for any $y\in\mathbb{R}^{n}$,

$\displaystyle(\operatorname{con}f_{\Lambda})(y)$
$\displaystyle=\inf\left\\{\sum_{i=1}^{n+1}\lambda_{i}f_{\Lambda}(y^{(i)}):\lambda_{i}\geq
0,\,\sum_{i=1}^{n+1}\lambda_{i}=1,\,y=\sum_{i=1}^{n+1}\lambda_{i}y^{(i)},\,y^{(i)}\in\mathbb{R}^{n}\right\\}$
$\displaystyle=\inf\left\\{\sum_{i=1}^{n+1}\lambda_{i}f(x^{(i)}):\lambda_{i}\geq
0,\,\sum_{i=1}^{n+1}\lambda_{i}=1,\,y=\sum_{i=1}^{n+1}\lambda_{i}\varphi(x^{(i)}),\,x^{(i)}\in\Omega\right\\}\geq\gamma.$

Let $x\in\Omega$. Since the only possible convex combination of the extreme
point $\varphi(x)$ from points $y^{(i)}\in\varphi(\Omega)$ is $\varphi(x)$
itself, we have

$\displaystyle(\operatorname{con}f_{\Lambda})(\varphi(x))$
$\displaystyle=\inf\left\\{\sum_{i=1}^{n+1}\lambda_{i}f_{\Lambda}(y^{(i)}):\lambda_{i}\geq
0,\,\sum_{i=1}^{n+1}\lambda_{i}=1,\,\varphi(x)=\sum_{i=1}^{n+1}\lambda_{i}y^{(i)},\,y^{(i)}\in\varphi(\Omega)\right\\}$
$\displaystyle=f_{\Lambda}(\varphi(x))=f(x).$

This shows that $\operatorname{con}f_{\Lambda}\circ\varphi=f$. Since
$\operatorname{con}f_{\Lambda}$ is bounded from below
$\operatorname{con}f_{\Lambda}$ is proper.

Let $f$ be lsc. Then $f_{\Lambda}$ inherits its lower semicontinuity from $f$:
Assume that $y\in\varphi(\Omega)$. Then there exists $x\in\Omega$ with
$y=\varphi(x)$ and we have due to the continuity of $\varphi$ and
$\varphi^{-1}$:

$\displaystyle\liminf_{y^{\prime}\to\varphi(x)}~{}f_{\Lambda}(y^{\prime})$
$\displaystyle:=\lim_{\delta\to
0^{+}}\inf~{}\\{f_{\Lambda}(y^{\prime}):y^{\prime}\in
B_{\delta}(\varphi(x))\\}$ $\displaystyle=\lim_{\delta\to
0^{+}}\inf~{}\\{f_{\Lambda}(\varphi(x^{\prime})):\varphi(x^{\prime})\in
B_{\delta}(\varphi(x))\\}=\lim_{\epsilon\to
0^{+}}\inf~{}\\{f(x^{\prime}):x^{\prime}\in B_{\epsilon}(x)\\}$
$\displaystyle=\liminf_{x^{\prime}\to x}~{}f(x^{\prime})\geq
f(x)=f_{\Lambda}(y).$

We invoke [33, Cor. 3.47] and deduce that $\operatorname{con}f_{\Lambda}$ is
proper, lsc and convex. In view of [33, Thm. 11.1] we have
$\operatorname{con}f_{\Lambda}=\operatorname{cl}\operatorname{con}f_{\Lambda}=f_{\Lambda}^{**}$.
∎

For a geometric intuition of this theorem we refer to Fig. 2.

For a piecewise polynomial discretization with degree at least $2$, due to
extremality, the primal discretized energy $\mathcal{F}_{\Lambda}$ restricted
to $\varphi(\Omega)^{\mathcal{V}}$ agrees with the original energy $F$. In
particular, this implies that an obtained Dirac solution
$(\varphi(x_{u}^{*}))_{u\in\mathcal{V}}$ of the discretization corresponds to
a solution of the original problem in the same way integer solutions of LP
relaxations are certificates of optimality for the corresponding ILP.

######  Proposition 3.12.

Let $\emptyset\neq\Omega\subset\mathbb{R}^{m}$. Let the metric $d$ be induced
by a norm and assume the space of linear functions on $\Omega$ is contained in
$\Lambda$. Furthermore, assume $\Lambda$ is spanned by an extremal curve
$\varphi:\Omega\to\mathbb{R}^{n}$ and $f_{u}$ lsc. Then for any
$y\in(\mathbb{R}^{n})^{\mathcal{V}}$, with $y_{u}=\varphi(x_{u})$,
$x_{u}\in\Omega$ the following identity holds true:

$\mathcal{F}_{\Lambda}(y)=F(x).$

In particular, whenever $y^{*}$ is a solution of problem (dR-P) such that
$y^{*}\in(\mathbb{R}^{n})^{\mathcal{V}}$, with $y_{u}^{*}=\varphi(x_{u}^{*})$,
for some $x_{u}^{*}\in\Omega$, $x^{*}$ is a solution of the original problem
(P).

###### Proof.

Let $d$ be induced by some norm $\lVert\cdot\rVert$ and denote its dual norm
by $\lVert\cdot\rVert_{*}$. As shown in Thm. 3.11, the unaries preserve the
original cost functions at $\varphi(\Omega)$. Hence it remains to show that
the for the pairwise costs it holds
$\sigma_{\mathcal{K}}((\nabla_{\Lambda}y)_{(u,v)})=\|x_{u}-x_{v}\|$. By
assumption $(\nabla_{\Lambda}y)_{(u,v)}=\varphi(x_{u})-\varphi(x_{v})$. We
rewrite

$\sigma_{\mathcal{K}_{\Lambda}}(\varphi(x_{u})-\varphi(x_{v}))=\sup_{\lambda\in\Lambda\cap\mathrm{Lip}_{d}(\Omega)}\lambda(x_{u})-\lambda(x_{v})\leq\|x_{u}-x_{v}\|.$

For any $x_{u},x_{v}\in\Omega$ we have

$\displaystyle\lVert
x_{u}-x_{v}\rVert=\sup_{\begin{subarray}{c}p\in\mathbb{R}^{m}:\|p\|_{*}\leq
1\end{subarray}}|\langle x_{u}-x_{v},p\rangle|=\langle
x_{u}-x_{v},p^{*}\rangle,$

where $p^{*}$ denotes the maximizer in the supremum, which exists due to the
compactness of the unit ball in a finite-dimensional space. Define the linear
function $\lambda^{*}:=\langle\cdot,p^{*}\rangle$ and note that by assumption
$\lambda^{*}\in\Lambda$. In addition we have shown that
$\lambda^{*}\in\mathrm{Lip}_{d}(\Omega)$ is $1$-Lipschitz. This implies
$\lVert
x_{u}-x_{v}\rVert=\lambda^{*}(x_{u})-\lambda^{*}(x_{v})\leq\sup_{\lambda\in\Lambda\cap\mathrm{Lip}_{d}(\Omega)}\lambda(x_{u})-\lambda(x_{v})\leq\lVert
x_{u}-x_{v}\rVert$ and hence equality holds. ∎

###  3.3 A Generalized Conjugacy Perspective

The previous results can be obtained from a generalized conjugacy point of
view. In particular, the convex conjugate of the lifted function
$f_{\Lambda}^{*}$ is comprised by the notion of $\Phi$-conjugacy, see [33, Ch.
11L*]:

######  Definition 3.13 (generalized conjugate functions).

Let $X$ and $Y$ be nonempty sets. Let $\Phi:X\times Y\to\overline{\mathbb{R}}$
be any function. Let $f:X\to\overline{\mathbb{R}}$. Then the $\Phi$-conjugate
of $f$ on $Y$ at $y\in Y$ is defined by

$\displaystyle f^{\Phi}(y):=\sup_{x\in X}\Phi(x,y)-f(x),$ (25)

and the $\Phi$-biconjugate of $f$ back on $X$ at $x\in X$ is given by

$\displaystyle f^{\Phi\Phi}(x):=\sup_{y\in Y}\Phi(x,y)-f^{\Phi}(y).$ (26)

We say that $f$ is a $\Phi$-envelope on $X$ if $f$ can be written in terms of
a pointwise supremum of a collection of elementary functions
$x\mapsto\Phi(x,y)-\alpha$, where $(\alpha,y)\in\overline{\mathbb{R}}\times Y$
is the parameter element.

Let $g:Y\to\overline{\mathbb{R}}$. Then, the $\Phi$-conjugate of $g$ on $X$ at
$x\in X$ is defined by

$\displaystyle g^{\Phi}(x):=\sup_{y\in Y}\Phi(x,y)-g(y),$ (27)

and the $\Phi$-biconjugate of $g$ back on $Y$ at $y\in Y$ is given by

$\displaystyle g^{\Phi\Phi}(y):=\sup_{x\in X}\Phi(x,y)-g^{\Phi}(x).$ (28)

We say that $g$ is a $\Phi$-envelope on $Y$ if $g$ can be written in terms of
a pointwise supremum of a collection of elementary functions
$y\mapsto\Phi(x,y)-\beta$, where $(\beta,x)\in\overline{\mathbb{R}}\times X$
is the parameter element.

For $X=\Omega$, $Y=\mathbb{R}^{n}$, and $\Phi(x,y)=\langle
y,\varphi(x)\rangle$, the convex conjugate $f_{\Lambda}^{*}$ is identical to
the $\Phi$-conjugate $f^{\Phi}$ of $f$, while its biconjugate
$f_{\Lambda}^{**}$ is the tightest lsc convex extension of $f_{\Lambda}$ to
$\operatorname{con}\varphi(\Omega)$. The $\Phi$-biconjugate $f^{\Phi\Phi}$ of
$f$ at a point $x\in\Omega$ is equal to the classical biconjugate of
$f_{\Lambda}$, evaluated at $\varphi(x)$, i.e.,
$f^{\Phi\Phi}=f_{\Lambda}^{**}\circ\varphi$ on $\Omega$, showing that the
$\Phi$-biconjugate is a convexely composite function and, therefore, it is
nonconvex in general. Actually, $\Phi$-conjugacy also comprises lifting to
measures via $\varphi(x)=\delta_{x}$, $\Phi$ the corresponding dual pairing
and $Y=\mathcal{C}(\Omega)$.

As a consequence of [33, Ex. 11.63], the considered $\Phi$-conjugacy can be
interpreted in terms of under-approximation by functions in $\Lambda$. In
analogy to the biconjugate $f^{**}$, which is the pointwise supremum of
affine-linear functions majorized by $f$, the $\Phi$-biconjugate
$f^{\Phi\Phi}$ is the pointwise supremum of functions in $\Lambda$ up to
constant translation majorized by $f$. This point of view also relates
$(f^{\Phi})^{*}$ and $f^{\Phi\Phi}$ by each other.

######  Remark 3.14.

The function $(f^{\Phi})^{*}=f_{\Lambda}^{**}$ is the pointwise supremum of
all affine-linear functions
$l_{\lambda,\beta}:=\langle\cdot,\lambda\rangle-\beta$ for which
$q_{\lambda,\beta}:=\Phi(\cdot,\lambda)-\beta$ is majorized by $f$. This can
be seen as follows: The Legendre–Fenchel conjugate can be characterized via
the identity
$(f^{\Phi})^{*}(y)=\sup_{(\lambda,\beta)\in\operatorname{epi}(f^{\Phi})}\langle
y,\lambda\rangle-\beta$. The observation now follows from the fact that
$(\lambda,\beta)\in\operatorname{epi}(f^{\Phi})$ if and only if
$q_{\lambda,\beta}$ is majorized by $f$.

The correspondence between $f_{\Lambda}^{**}$ and $(f^{\Phi})^{*}$ and between
the minorizers $l_{\lambda,\beta}$ and $q_{\lambda,\beta}$ is illustrated in
Fig. 2. Note that this is closely related to the idea of feature maps
$\varphi$ in linear classifiers.

Thm. 3.11 identifies all lsc functions $f$ as $\Phi$-envelopes whenever
$\varphi$ is extremal:

######  Corollary 3.15.

Let $\Omega\subset\mathbb{R}^{m}$ be nonempty and compact and let
$f:\Omega\to\mathbb{R}$ be bounded from below. Furthermore, let
$\varphi:\Omega\to\mathbb{R}^{n}$ be an extremal curve. Then we have

$\displaystyle f^{\Phi\Phi}=f_{\Lambda}^{**}\circ\varphi=\operatorname{cl}f,$
(29)

on $\Omega$.

###### Proof.

Since $f$ is finite-valued and bounded from below on $\Omega$ we have
$(\operatorname{cl}f)(x)>-\infty$ for $x\in\Omega$ and therefore
$\operatorname{cl}f$ is finite-valued. By [33, Ex. 11.63] $f^{\Phi\Phi}$ is
the largest $\Phi$-envelope below $f$. Since $\varphi$ is continuous relative
to $\Omega$, $f^{\Phi\Phi}=f_{\Lambda}^{**}\circ\varphi$ is lsc relative to
$\Omega$. Since $\operatorname{cl}f$ is the largest lsc function below $f$ we
have $f^{\Phi\Phi}\leq\operatorname{cl}f$. Since $\operatorname{cl}f\leq f$ we
also have $(\operatorname{cl}f)^{\Phi\Phi}\leq f^{\Phi\Phi}$. Invoking Thm.
3.11 we have

$\operatorname{cl}f=(\operatorname{cl}f)^{\Phi\Phi}\leq
f^{\Phi\Phi}\leq\operatorname{cl}f,$

on $\Omega$. Therefore $f^{\Phi\Phi}=\operatorname{cl}f$ on $\Omega$. ∎

Up to the presence of the compact set $\Omega$, this result generalizes the
basic quadratic transform [33, Ex. 11.66] originally due to [27, Prop. 3.4]
(for lsc functions only), which is obtained by choosing
$\varphi(x)=(x_{1},x_{2},\ldots,x_{m},\|x\|^{2})$.

##  4 A Tractable Conic Program for MRFs

###  4.1 Nonnegativity and Moments

After discretization a next step to obtain a practical implementation is to
derive finite characterizations of the lifted biconjugates $f_{\Lambda}^{**}$
and the constraint set $\mathcal{K}_{\Lambda}$. We will show that the
formulations can be rewritten in terms of a semi-infinite conic program which
can be implemented using semidefinite programming in the piecewise polynomial
case.

For now let $f_{u}\in\Lambda$. Eq. (21) then shows, that the challenging part
is to characterize the moment space $\mathcal{P}_{\Lambda}$: The following
result shows that up to normalization, $\mathcal{P}_{\Lambda}$ can be written
in terms of the dual cone of the cone of functions in $\Lambda$ that are
nonnegative on $\Omega$.

######  Lemma 4.1.

Let $(\mathcal{M}_{\Lambda})_{+}$ be the cone of moments of nonnegative
measures as defined in (12) and let $\mathcal{N}_{\Lambda}$ be the cone of the
coefficients of the functions in
$\Lambda=\langle\varphi_{0},\ldots,\varphi_{n}\rangle$ that are nonnegative on
$\Omega$ defined as:

$\displaystyle\mathcal{N}_{\Lambda}:=\\{p\in\mathbb{R}^{n+1}:\langle
p,\varphi(x)\rangle\geq 0,\;\forall\,x\in\Omega\\}.$ (30)

Then $(\mathcal{M}_{\Lambda})_{+}$ is equal to $\mathcal{N}_{\Lambda}^{*}$,
where $\mathcal{N}_{\Lambda}^{*}$ denotes the dual cone of
$\mathcal{N}_{\Lambda}$.

If, in addition, $\varphi_{0}\equiv 1$ we also have
$\mathcal{P}_{\Lambda}=\\{y\in(\mathcal{M}_{\Lambda})_{+}\colon y_{0}=1\\}$.

###### Proof.

Let $y\in(\mathcal{M}_{\Lambda})_{+}$. This means there exists
$\mu\in\mathcal{M}_{+}(\Omega)$ such that
$y_{k}:=\int_{\Omega}\varphi_{k}(x)\,\mathrm{d}\mu(x)$. Let
$p\in\mathcal{N}_{\Lambda}$. Because of $\langle
p,\varphi(\cdot)\rangle\in\Lambda\subset\mathcal{C}(\Omega)$ and $\langle
p,\varphi(x)\rangle\geq 0$ for all $x\in\Omega$ and
$\mu\in\mathcal{M}_{+}(\Omega)$ is a nonnegative measure it holds $\langle
p,y\rangle=\sum_{k=0}^{n}p_{k}\int_{\Omega}\varphi_{k}(x)\,\mathrm{d}\mu(x)=\int_{\Omega}\langle
p,\varphi(x)\rangle\,\mathrm{d}\mu(x)\geq 0$ for all $x\in\Omega$. Since
$p\in\mathcal{N}_{\Lambda}$ was an arbitrary choice from
$\mathcal{N}_{\Lambda}$ we have $y\in\mathcal{N}_{\Lambda}^{*}$.

Next we show $(\mathcal{M}_{\Lambda})_{+}^{*}\subseteq\mathcal{N}_{\Lambda}$
as this implies
$\mathcal{N}_{\Lambda}^{*}\subseteq(\mathcal{M}_{\Lambda})_{+}^{**}=(\mathcal{M}_{\Lambda})_{+}$,
where the last equality holds since $(\mathcal{M}_{\Lambda})_{+}$ is closed
and convex. Take $p\in(\mathcal{M}_{\Lambda})_{+}^{*}$. Let $x\in\Omega$. Now
we choose $y\in(\mathcal{M}_{\Lambda})_{+}$ such that
$y_{k}=\langle\delta_{x},\varphi_{k}\rangle=\varphi_{k}(x)$. Then
$p\in(\mathcal{M}_{\Lambda})_{+}^{*}$ implies that $\langle p,y\rangle\geq 0$.
Since the choice $x\in\Omega$ was arbitrary we have $\langle
p,\varphi(x)\rangle\geq 0$ for all $x\in\Omega$ and therefore
$p\in\mathcal{N}_{\Lambda}$.

Finally, $\mathcal{P}_{\Lambda}=\\{y\in(\mathcal{M}_{\Lambda})_{+}\colon
y_{0}=1\\}$ follows from the fact that $\mu\in\mathcal{M}(\Omega)$ is an
element of $\mathcal{P}(\Omega)$ if and only if
$\mu\in\mathcal{M}_{+}(\Omega)$ and $\langle\mu,\varphi_{0}\rangle=1$. ∎

Before we specialize $\Lambda$ to the space of polynomials, we derive a cone
programming formulation of the Lipschitz constraints
$\lambda_{(u,v)}\in\mathrm{Lip}_{d}(\Omega)$. Here we restrict $\Omega=[a,b]$
to be a compact interval. In the following, we provide an implementation for
two specific metrics. As before, this boils down to nonnegativity of
functions: Firstly, we consider total variation regularization, i.e.
$d(x,y)=|x-y|$: Assume that $\Lambda$ is closed under differentiation, i.e.,
$\varphi$ is differentiable and $\varphi_{k}^{\prime}\in\Lambda$. Then, the
condition $\lambda_{(u,v)}\in\mathrm{Lip}_{d}([a,b])$ can be phrased in terms
of the constraints $-1\leq\lambda_{(u,v)}^{\prime}(x)\leq 1$ for all
$x\in[a,b]$, where $\lambda_{(u,v)}^{\prime}$ is the derivative of
$\lambda_{(u,v)}$. Equivalently, this means that the coefficients of the
functions $1+\lambda_{(u,v)}^{\prime}$ and $1-\lambda_{(u,v)}^{\prime}$ are in
$\mathcal{N}_{\Lambda}$. Secondly, we consider Potts regularization, i.e.,
$d(x,y)=\llbracket x=y\rrbracket$. Since $\lambda_{(u,v)}$ is a univariate
function, and constant terms in the dual variable do not matter, the condition
$\lambda_{(u,v)}\in\mathrm{Lip}_{d}([a,b])$ can be equivalently phrased as
$0\leq\lambda_{(u,v)}\leq 1$, see, [43, Ex. 1.17]. Equivalently, this means
that the coefficients of the functions $\lambda_{(u,v)}$ and
$1-\lambda_{(u,v)}$ are in $\mathcal{N}_{\Lambda}$.

###  4.2 Semidefinite Programming and Nonnegative Polynomials

As we have seen in the previous section, an important ingredient for a
tractable formulation is the efficient characterization of nonnegativity of
functions in a finite-dimensional subspace
$\Lambda\subset\mathcal{C}(\Omega)$. A promising choice of $\Lambda$ in that
regards is the space of polynomials. Indeed, the characterization of
nonnegativity of polynomials is a fundamental problem in _convex algebraic
geometry_ surveyed in [5]: Let $\mathbb{R}[x_{1},\ldots,x_{m}]$ denote the
ring of possibly multivariate polynomials with
$p\in\mathbb{R}[x_{1},\ldots,x_{m}]$. Then $p=\sum_{\alpha\in
I}p_{\alpha}x^{\alpha}$ for monomials $x^{\alpha}$. Let $\deg p$ denote its
degree. A key result from real algebraic geometry is the Positivstellensatz
due to [17] and [41] refined in [39] and [29]. It characterizes polynomials
$p\in\mathbb{R}[x_{1},\ldots,x_{m}]$ that are positive on semi-algebraic sets
$\Omega$, i.e. $p(x)>0$ for all $x\in\Omega$, where $\Omega$ is defined in
terms of polynomial inequalities.

Key to such results is a _certificate_ of nonnegativity of the polynomial $p$
that involves _sum-of-squares_ (SOS) multipliers $q$, where $q$ is SOS if
$q(x)=\sum_{i=1}^{N}q_{i}^{2}(x)$ for polynomials
$q_{i}\in\mathbb{R}[x_{1},\ldots,x_{m}]$. For intervals $X=[a,b]$, thanks to
[5, Thm. 3.72] originally due to [28, Cor. 2.3] we have following result:

######  Lemma 4.2.

Let $a<b$. Then the univariate polynomial $p\in\mathbb{R}[x]$ is nonnegative
on $[a,b]$ if and only if it can be written as

$\displaystyle p(x)=\begin{cases}s(x)+(x-a)\cdot(b-x)\cdot t(x)&\text{if $\deg
p$ is even,}\\\ (x-a)\cdot s(x)+(b-x)\cdot t(x)&\text{if $\deg p$ is
odd,}\end{cases}$ (31)

where $s,t\in\mathbb{R}[x]$ are sum of squares. If $\deg p=2n$, then we have
$\deg s\leq 2n$, $\deg t\leq 2n-2$, while if $\deg p=2n+1$, then $\deg s\leq
2n,\deg t\leq 2n$.

Remarkably, the above result provides us with explicit upper bounds of the
degrees of the SOS multipliers $s$ and $t$ that are important to derive a
practical implementation: Then the SOS constraints can be formulated in terms
of semidefinite and affine inequalities: We adopt [5, Lem. 3.33] and [5, Lem.
3.34]:

######  Lemma 4.3.

A univariate polynomial $p\in\mathbb{R}[x]$ with $\deg p=2n$, $n\geq 0$ is SOS
if and only if there exists a positive semidefinite matrix
$Q\in\mathbb{R}^{n+1\times n+1}$ such that

$\displaystyle p_{k}=\sum_{\begin{subarray}{c}0\leq i,j\leq n,\\\
i+j=k\end{subarray}}Q_{ij},\quad\forall\,0\leq k\leq 2n.$ (32)

Invoking the results above SDP-duality yields the following compact
representation of $\mathcal{M}_{\Lambda}$:

######  Lemma 4.4.

Let $n\geq 0$. For odd degree $2n+1$, $y\in\mathcal{M}_{\Lambda}$ if and only
if

$\displaystyle bM_{0,n}(y)\succeq M_{1,n}(y)\succeq aM_{0,n}(y),$ (33)

for Hankel matrices

$\displaystyle M_{i,n}(y):=\begin{bmatrix}y_{i}&y_{i+1}&\ldots&y_{i+n}\\\
y_{i+1}&y_{i+2}&\ldots&y_{i+n+1}\\\ \vdots&&\ldots&\vdots\\\
y_{i+n}&&\ldots&y_{i+2n}\end{bmatrix}.$ (34)

For even degree $2n$, $y\in\mathcal{M}_{\Lambda}$ if and only if

$\displaystyle M_{0,n}(y)\succeq 0,$ (35)
$\displaystyle\left(a+b\right)M_{1,n-1}(y)-abM_{0,n-1}(y)\succeq
M_{2,n-1}(y).$ (36)

###### Proof.

Follows by Lem. 4.2 and Lem. 4.3 invoking elementary SDP-duality. ∎

###  4.3 A First-Order Primal-Dual Algorithm

We are now ready to describe the algorithm for solving the resulting
semidefinite program. We first consider the case $f_{u}\in\Lambda$ and
$\Lambda$ is the space of univariate polynomials. We propose to use the PDHG
[9] algorithm, as it can exploit the partially separable structure of our SDP.
The primal-dual algorithm optimizes the problem (dR-P) via alternating
projected gradient descent/ascent steps applied to the saddle-point
formulation of (dR-P):

$\displaystyle\min_{y\in(\mathcal{P}_{\Lambda})^{\mathcal{V}}}~{}\max_{p\in(\mathcal{K}_{\Lambda})^{\mathcal{E}}}~{}~{}\langle
y,f-\operatorname{Div}_{\Lambda}p\rangle,$ (37)

which is obtained by expanding the support function in Problem (dR-P) and
substituting the expression (21) for the lifted biconjugates. In each
iteration the algorithm performs a projected gradient ascent step in the dual
$p$ followed by a projected gradient descent step in the primal variable $y$.
Subsequently it performs an extrapolation step in the primal. The projections
onto the sets $(\mathcal{P}_{\Lambda})^{\mathcal{V}}$ and
$(\mathcal{K}_{\Lambda})^{\mathcal{E}}$ are separable and can therefore be
carried out in parallel on a GPU using the SDP characterizations derived
above. For practicality, we introduce additional auxiliary variables and
linear constraints to decouple the affine constaints (32) and the SDP
constraints. The projection operator of the semidefinite cone can then be
solved using an eigenvalue decomposition.

###  4.4 Piecewise Polynomial Duals and Nonlinear Lifted Biconjugates

The polynomial discretization can be extended by means of a continuous
piecewise polynomial representation of the dual variables resulting in a
possibly more accurate approximation of the dual subspace $\Lambda$. Then,
both, nonnegativity and Lipschitz continuity can be enforced on each piece
$\Omega_{k}$ individually. Continuity of the piecewise polynomial dual
variables can be enforced via linear constraints. The corresponding primal
variable $y$ belongs to
$y\in(\mathcal{M}_{\Lambda})_{+}\times(\mathcal{M}_{\Lambda})_{+}\times\cdots\times(\mathcal{M}_{\Lambda})_{+}$.
Then the restriction that $y$ is a moment vector of a probability measure
supported on the whole space $\Omega$ yields an additional sum-to-one
constraint on the $0^{\text{th}}$ moments $1=\sum_{k=1}^{K}y_{k,0}$.

Another issue to address is when $f_{u}\not\in\Lambda$ which results in a
nonlinear lifted biconjugate over the moment-space as in Fig. 2.

The formulation which is derived next addresses both: In particular it allows
one to choose $\Lambda$ independently from $f_{u}$ which can even be
discontinuous, as long as $f_{u}$ has a piecewise polynomial structure. Key to
the formulation is to rewrite the inner minimum in the dual formulation dR-D
exploiting a duality between nonnegativity and minimization of functions:

$\displaystyle 1$$\displaystyle 2$$\displaystyle 3$$\displaystyle 4$Number of
pieces $\displaystyle
k$$\displaystyle-240$$\displaystyle-220$$\displaystyle-200$$\displaystyle-180$$\displaystyle-160$$\displaystyle-140$EnergyTV
regularizationdeg$\displaystyle=1$deg$\displaystyle=2$deg$\displaystyle=3$deg$\displaystyle=4$

(a)

$\displaystyle 1$$\displaystyle 2$$\displaystyle 3$$\displaystyle 4$Number of
pieces $\displaystyle
k$$\displaystyle-250$$\displaystyle-200$$\displaystyle-150$$\displaystyle-100$$\displaystyle-50$$\displaystyle
0$EnergyTV
regularizationdeg$\displaystyle=1$deg$\displaystyle=2$deg$\displaystyle=3$deg$\displaystyle=4$

(b)

$\displaystyle 1$$\displaystyle 2$$\displaystyle 3$$\displaystyle 4$Number of
pieces $\displaystyle k$$\displaystyle 10^{1}$$\displaystyle 10^{2}$Primal-
Dual gapTV
regularizationdeg$\displaystyle=1$deg$\displaystyle=2$deg$\displaystyle=3$deg$\displaystyle=4$

(c)

$\displaystyle 1$$\displaystyle 2$$\displaystyle 3$$\displaystyle 4$Number of
pieces $\displaystyle
k$$\displaystyle-250$$\displaystyle-225$$\displaystyle-200$$\displaystyle-175$$\displaystyle-150$$\displaystyle-125$$\displaystyle-100$$\displaystyle-75$EnergyPotts
regularizationdeg$\displaystyle=1$deg$\displaystyle=2$deg$\displaystyle=3$deg$\displaystyle=4$

(d)

Figure 3: Primal and dual energies for MAP-inference in a continuous MRF with
TV regularization using a piecewise polynomial hierarchy of dual variables.
(a) shows the dual energy for TV. In (b) the dashed lines correspond to the
primal energy at the rounded solution and the solid lines correspond to the
dual energy. (c) shows the gap between the nonconvex primal energy at the
rounded solution and the dual energy for TV regularization. (d) shows the dual
energies for Potts regularization.

Let $\Omega=[a,b]$, $a<b$ be a compact interval. Let
$a=t_{1}<t_{2}<t_{3},\ldots,<t_{K+1}=b$ be a sequence of knots, where
$\Omega_{k}:=[t_{k},t_{k+1}]$. Let $\Theta$ be the space of univariate
polynomials with some maximum degree $n$. Let $f:\Omega\to\mathbb{R}$ be a
possibly discontinuous lsc piecewise polynomial function defined by
$f(x)=\min_{1\leq k\leq K}f_{k}(x)+\iota_{\Omega_{k}}(x)$ with
$f_{k}\in\Theta$, i.e., $f_{k}(x)=\langle\varphi(x),w_{k}\rangle$ for
coefficients $w_{k}\in\mathbb{R}^{n+1}$, where $\varphi_{0}\equiv 1$. First
observe the following duality between nonnegativity and minimization of a lsc
function:

$\min_{x\in\Omega}~{}f(x)=\max_{q\in\mathbb{R}}~{}q-\iota_{\mathcal{N}}(f-q),$

where we denote by
$\mathcal{N}(\Omega)=\\{\lambda:\Omega\to\mathbb{R}:\lambda(x)\geq
0,\forall\;x\in\Omega\\}$ the cone of nonnegative lsc functions on $\Omega$.
Then we obtain for $Aq=(e_{0}q,\ldots,e_{0}q)$, where
$e_{0}=(1,0,\ldots,0)\in\mathbb{R}^{n+1}$ is the $0^{\text{th}}$ unit vector:

$\displaystyle\min_{x\in\Omega}~{}f(x)$
$\displaystyle=\max_{q\in\mathbb{R}}~{}q-\iota_{\mathcal{N}(\Omega)}(f-q)$
(38)
$\displaystyle=\max_{q\in\mathbb{R}}~{}q-\sum_{k=1}^{K}\iota_{\mathcal{N}_{\Theta}}(w_{k}-A_{k}q).$
(39)

Fenchel–Rockafellar duality then yields:

$\displaystyle\min_{x\in\Omega}~{}f(x)$
$\displaystyle=\min_{y\in(\mathbb{R}^{n+1})^{K}}\iota_{\\{1\\}}(A^{*}y)+\sum_{k=1}^{K}\langle
y_{k},w_{k}\rangle+\iota_{(\mathcal{M}_{\Theta})_{+}}(y_{k})$
$\displaystyle=\min_{\begin{subarray}{c}y\in((\mathcal{M}_{\Theta})_{+})^{K}\\\
\sum_{k=1}^{K}y_{k,0}=1\end{subarray}}\sum_{k=1}^{K}\langle
y_{k},w_{k}\rangle.$

This formulation can be substituted in the dual problem (dR-D) and we obtain

$\displaystyle\sup_{p\in((\mathbb{R}^{n})^{K})^{\mathcal{E}}}\sum_{u\in\mathcal{V}}~{}\min_{\begin{subarray}{c}y\in((\mathcal{M}_{\Theta})_{+})^{K}\\\
\sum_{k=1}^{K}y_{k,0}=1\end{subarray}}\sum_{k=1}^{K}\langle
y_{k},(w-\operatorname{Div}_{\Lambda}p)_{u,k}\rangle-\sum_{e\in\mathcal{E}}\iota_{\mathcal{K}}(p_{e}).$
(40)

Here the dual variables $\lambda$ are chosen such that
$(\operatorname{Div}_{\Lambda}p)_{u}$ represents a piecewise polynomial with
knots $a=t_{1}<t_{2}<t_{3}<\cdots<t_{K+1}=b$ such that for each piece we have
$\langle(\operatorname{Div}_{\Lambda}p)_{u,k},\varphi(\cdot)\rangle\in\Theta$.
Note that this does not require $\Lambda$ to be equal to the whole space of
continuous piecewise polynomials of degree $n$. Indeed, $\Lambda$ can be a
subspace thereof which covers the case where $f_{u}\not\in\Lambda$.

##  5 Numerical Experiments

###  5.1 Empirical Convergence Study

Left image stereo pair | Standard $k=30$ | $k=5,\deg=1$ | $k=5,\deg=7$
---|---|---|---
|  |  |
| rounded $24611.51$ | rounded $22283.25$ | rounded $19428.49$
|  | dual $16227.80$ | dual $17472.13$
|  |  |
| rounded $17510.55$ | rounded $15027.55$ | rounded $13380.71$
|  | dual $10962.53$ | dual $12008.45$
Figure 4: Stereo disparity estimation from a stereo image pair: Left: Standard
MRF/OT discretization implemented using a continuous piecewise linear under-
approximation for the unaries and piecewise linear duals. middle: piecewise
linear duals. right: piecewise polynomial duals. The visual appearence of the
solution to the standard MRF/OT discretization shows a strong grid bias. The
dual energy gap increases for increasing the degree and/or the number of
pieces. Likewise the energy at the rounded solution decreases.

In this first experiment we evaluate the local marginal polytope relaxation of
the MRF formulation (P) using a piecewise polynomial hierarchy of dual
variables. We choose a vertex grid of size $16\times 16$. We fix a random
polynomial data term of degree $4$ at each vertex by fitting a random sample
of data points. To obtain a high-accuracy solution we solve the primal SDP
formulation corresponding to the saddle-point formulation (40) with
MOSEK111https://www.mosek.com/products/academic-licenses. For recovering a
primal solution at each vertex $u$ we compute the mode w.r.t. the
$0^{\text{th}}$ moments to select the best interval denoted by
$k^{*}=\operatorname*{arg\,max}_{1\leq k\leq K}(y_{u})_{k,0}$. Then we compute
the mean of the discretized measure corresponding to the $(k^{*})^{\text{th}}$
interval as $x_{u}=(y_{u})_{k^{*},1}$.

Fig. 3 visualizes the primal and dual energies for varying degrees and/or
number of pieces of the dual variable. While the dual energy strictly
increases with higher degrees and/or number of pieces the primal energy is
evaluated at the rounded solution and therefore does not strictly decrease in
general. While for TV increasing the degree vs. increasing the number of
pieces (for $K\cdot\deg$ constant) leads to similar performance, for Potts, in
many situations, increasing the degree leads to larger dual energies, e.g.,
consider $\deg=4,K=1$ vs. $\deg=1,K=4$, red curve vs. blue curve in Fig. 3(d).
In further experiments, we observed, that this holds in particular when the
structure of the dual variables and the unaries match, i.e.,
$f_{u},\lambda_{e}\in\Lambda$. Note that for Potts, since the dual variables
are uniformly bounded on $\Omega$ and the derivative can be unbounded we drop
the continuity constraint which leads to a more compact formulation and larger
dual energies.

Dual energies
---
$\deg$ | $K=1$ | $K=3$ | $K=5$
$1$ | 14180.08 | 15733.71 | 16227.80
$2$ | 15052.18 | 16430.97 | 16773.75
$3$ | 15601.89 | 16778.87 | 17055.25
$4$ | 15938.92 | 16998.63 | 17235.21
$5$ | 16191.40 | 17147.92 | 17346.49
$6$ | 16369.95 | 17243.47 | 17422.26
$7$ | 16480.22 | 17308.91 | 17472.13

Energies rounded
---
$\deg$ | $K=1$ | $K=3$ | $K=5$
$1$ | 28982.45 | 25005.05 | 22283.26
$2$ | 31038.14 | 22525.12 | 21049.32
$3$ | 28505.41 | 21500.65 | 20428.34
$4$ | 27255.48 | 20841.62 | 20049.04
$5$ | 25795.56 | 20344.17 | 19764.23
$6$ | 24081.51 | 20032.77 | 19559.00
$7$ | 23142.67 | 19869.71 | 19428.50

Table 1: Energies for stereo matching Motorcycle. Left: Dual energies. Right:
Primal energies at the rounded solution.

###  5.2 Stereo Matching

In this experiment we consider stereo matching using the anisotropic
relaxation (R-P). We consider the Motorcycle and the Piano image pairs from
the Middlebury stereo benchmark [38]. We downsample the images by factor $4$.
The disparity cost term is first calculated using 135 discrete disparities
obtained by shifting the images by the corresponding amount of pixels and
comparing the image gradients. Then, the cost dataterm is approximated from
below in terms of a continuous piecewise cubic polynomial $f_{u}$ using 30
pieces at each $u\in\mathcal{V}$. In Fig. 4 we compare the standard MRF/OT
discretization as described in Ex. 3.3 with our framework using piecewise
linear and piecewise quadratic dual variables. The standard MRF/OT
discretization is equivalent to a piecewise linear approximation of the data
term with piecewise linear duals in our framework. For a fair comparison we
use a continuous piecewise linear under-approximation of the continuous
piecewise cubic polynomial under-approximation $f_{u}$. As the resulting
optimization problem is large-scale we solve the saddle-point formulation (37)
with PDHG [9] as described in Sec. 4.3 using the GPU-based PDHG framework
prost222https://github.com/tum-vision/prost. In contrast to the previous
experiment which uses a combined mode and mean rounding procedure we found the
plain mean of the discretized measure to produce better results on real data:
More explicitly we recover a solution according to
$x_{u}=\sum_{k=1}^{K}t_{k}(y_{u})_{k,0}$ at each vertex $u\in\mathcal{V}$. In
Tab. 1 we compare both, dual and nonconvex primal energies, for a larger
hierarchy of dual subspaces.

##  6 Discussion

We presented a method to reduce duality gaps in the Lagrangian relaxation of
the MAP-inference problem in a continuous MRF taking a nonlinear optimization-
driven approach. Our theoretical contribution identifies _extremality_ of the
lifting as a key component, as it leads to convexifications which do not
discard any information for a wide range of nonconvex cost functions. Using
results from convex algebraic geometry we provide, to our knowledge, the first
tractable formulation of the polynomial discretization in terms of
semidefinite programming. We have provided a parallel implementation of a
first-order primal dual algorithm on a GPU which can handle large problems.
Indeed, the approach of [18, 19] applied directly to the original problem (P)
(with $f_{u},f_{uv}$ polynomial) attempts to solve the _full_ marginal
polytope relaxation which is tight but intractable for large $\mathcal{V}$ as
the number of coupling moments grows exponentially with $|\mathcal{V}|$. In
contrast, our framework applies to the _local_ marginal polytope relaxation
which is not tight in general but leads to a tractable formulation as it
exploits the sparse structure of the optimization problem.

As revealed by our experiments, both increasing the number of pieces and the
degree of the dual variables successfully reduces the nonconvex duality gap.
In particular, the algorithm is applied to the stereo matching problem between
two images, showing significant improvements over piecewise constant or
piecewise linear discretizations in the dual which demand a higher number of
samples.

For total variation regularized problems, our theory suggests that the duality
gap vanishes like $\mathcal{O}(1/(K\cdot\sqrt{\deg}))$. The faster convergence
in the number of pieces $K$ is also confirmed by our experiments. However,
setting the degree $\deg\geq 2$ has the attractive theoretical property that
all information of the original nonconvex cost is preserved. For Potts
regularization, our experiments indicate that increasing the degree leads to
better results than increasing the pieces.

Going beyond the piecewise-linear setting might be particularly promising for
vector- or manifold-valued label-spaces with $\dim(\Omega)>1$. Labeling
problems in such higher-dimensonal settings were recently considered in [21,
24, 46]. There, due to the piecewise-linear dual variables, the complexity
grows exponentially with $\dim(\Omega)$. Moreover, as discussed in [46], the
piecewise-approach inherently requires an approximation of $\Omega$ by a
triangulated manifold. Building upon the techniques developed in this paper,
an interesting direction for future work is to investigate whether these two
drawbacks can be overcome by considering suitably defined approximation spaces
for the dual variables in vector- or manifold-valued settings.

## Acknowledgments

We would like to thank Johannes Milz, Peter Ochs and Jan-Hendrik Lange for
their valuable feedback on an early version of this manuscript.

## References

  * [1] E. J. Balder.  An extension of duality-stability relations to nonconvex optimization problems.  SIAM J. Control Optim., 15(2):329–343, 1977.
  * [2] R. Bergmann, R. H. Chan, R. Hielscher, J. Persch, and G. Steidl.  Restoration of manifold-valued images by half-quadratic minimization.  arXiv preprint arXiv:1505.07029, 2015.
  * [3] F. Bernard, F. R. Schmidt, J. Thunberg, and D. Cremers.  A combinatorial solution to non-rigid 3d shape-to-image matching.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017.
  * [4] D. Bertsekas.  Nonlinear Programming.  Athena Scientific, 1999.
  * [5] G. Blekherman, P. A. Parrilo, and R. R. Thomas.  Semidefinite Optimization and Convex Algebraic Geometry.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2012\.
  * [6] B. Brown, D. Elliott, and D. Paget.  Lipschitz constants for the Bernstein polynomials of a Lipschitz continuous function.  Journal of approximation theory, 49(2):196–199, 1987.
  * [7] C. Carathéodory.  Über den Variabilitätsbereich der Fourier’schen Konstanten von positiven harmonischen Funktionen.  Rendiconti Del Circolo Matematico di Palermo (1884-1940), 32(1):193–217, 1911.
  * [8] N. L. Carothers.  A short course on approximation theory.  1998\.
  * [9] A. Chambolle and T. Pock.  A first-order primal-dual algorithm for convex problems with applications to imaging.  Journal of mathematical imaging and vision, 40(1):120–145, 2011\.
  * [10] Q. Chen and V. Koltun.  Full flow: Optical flow estimation by global optimization over regular grids.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2016.
  * [11] T. Chen, J.-B. Lasserre, V. Magron, and E. Pauwels.  Semialgebraic optimization for Lipschitz constants of ReLU networks.  In Conference on Neural Information Processing Systems, 2020.
  * [12] D. J. Crandall, A. Owens, N. Snavely, and D. P. Huttenlocher.  SfM with MRFs: Discrete-continuous optimization for large-scale structure from motion.  IEEE Trans. Pattern Anal. Mach. Intell. (PAMI), 35(12):2841–2853, 2012.
  * [13] C. Domokos, F. R. Schmidt, and D. Cremers.  MRF optimization with separable convex prior on partially ordered labels.  In European Conference on Computer Vision (ECCV), 2018.
  * [14] A. Fix and S. Agarwal.  Duality and the continuous graphical model.  In European Conference on Computer Vision (ECCV), 2014.
  * [15] L. V. Kantorovich.  Mathematical methods of organizing and planning production.  Management Science, 6(4):366–422, 1960.
  * [16] J. Kappes, B. Andres, F. Hamprecht, C. Schnorr, S. Nowozin, D. Batra, S. Kim, B. Kausler, J. Lellmann, N. Komodakis, et al.  A comparative study of modern inference techniques for discrete energy minimization problems.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2013.
  * [17] J.-L. Krivine.  Anneaux préordonnés.  Journal d’analyse mathématique, 12(1):307–326, 1964.
  * [18] J. B. Lasserre.  Global optimization with polynomials and the problem of moments.  SIAM Journal on Optimization, 11(3):796–817, 2001.
  * [19] J. B. Lasserre.  Semidefinite programming vs. LP relaxations for polynomial programming.  Mathematics of Operations Research, 27(2):347–360, 2002.
  * [20] F. Latorre, P. T. Y. Rolland, and V. Cevher.  Lipschitz constant estimation for neural networks via sparse polynomial optimization.  In 8th International Conference on Learning Representations, 2020\.
  * [21] E. Laude, T. Möllenhoff, M. Moeller, J. Lellmann, and D. Cremers.  Sublabel-accurate convex relaxation of vectorial multilabel energies.  In European Conference on Computer Vision (ECCV), 2016.
  * [22] C. Lemaréchal.  Lagrangian relaxation.  In Computational Combinatorial Optimization, pages 112–156. Springer, 2001.
  * [23] T. Möllenhoff and D. Cremers.  Sublabel-accurate discretization of nonconvex free-discontinuity problems.  In International Conference on Computer Vision (ICCV), 2017.
  * [24] T. Möllenhoff and D. Cremers.  Lifting vectorial variational problems: a natural formulation based on geometric measure theory and discrete exterior calculus.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2019.
  * [25] J. J. Moreau.  Inf-convolution, sous-additivité, convexité des fonctions numériques.  Journal de Mathématiques Pures et Appliquées, (47):33–41, 1970.
  * [26] J. Peng, T. Hazan, D. McAllester, and R. Urtasun.  Convex max-product algorithms for continuous MRFs with applications to protein folding.  In International Conference on Machine Learning (ICML), 2011.
  * [27] R. Poliquin.  Subgradient monotonicity and convex functions.  Nonlinear Analysis: Theory, Methods & Applications, 14(4):305–317, 1990.
  * [28] V. Powers and B. Reznick.  Polynomials that are positive on an interval.  Transactions of the American Mathematical Society, 352(10):4677–4692, 2000.
  * [29] M. Putinar.  Positive polynomials on compact semi-algebraic sets.  Indiana University Mathematics Journal, 42(3):969–984, 1993.
  * [30] R. T. Rockafellar.  Duality and stability in extremum problems involving convex functions.  Pacific J. Math., 21(1):167–187, 1967.
  * [31] R. T. Rockafellar.  Convex Analysis.  Princeton University Press, New Jersey, 1970.
  * [32] R. T. Rockafellar.  Augmented Lagrange multiplier functions and duality in nonconvex programming.  SIAM J. Control Optim., 12(2):268–285, 1974.
  * [33] R. T. Rockafellar and R. J.-B. Wets.  Variational Analysis.  Springer, New York, 1998.
  * [34] W. Rudin.  Functional analysis.  Internat. Ser. Pure Appl. Math, 1991.
  * [35] N. Ruozzi.  Exactness of approximate MAP inference in continuous MRFs.  In Advances in Neural Information Processing Systems (NeurIPS), 2015\.
  * [36] M. Salzmann.  Continuous inference in graphical models with polynomial energies.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2013.
  * [37] F. Santambrogio.  Optimal Transport for Applied Mathematicians.  Birkhäuser, New York, 2015.
  * [38] D. Scharstein, H. Hirschmüller, Y. Kitajima, G. Krathwohl, N. Nešić, X. Wang, and P. Westling.  High-resolution stereo datasets with subpixel-accurate ground truth.  In German conference on pattern recognition, pages 31–42. Springer, 2014.
  * [39] K. Schmüdgen.  The $k$-moment problem for compact semi-algebraic sets.  Mathematische Annalen, 289(1):203–206, 1991.
  * [40] M. Simões, A. Themelis, and P. Patrinos.  Lasry-lions envelopes and nonconvex optimization: A homotopy approach.  arXiv preprint arXiv:2103.08533, 2021.
  * [41] G. Stengle.  A Nullstellensatz and a Positivstellensatz in semialgebraic geometry.  Mathematische Annalen, 207(2):87–97, 1974.
  * [42] H. Trinh and D. McAllester.  Particle-based belief propagation for structure from motion and dense stereo vision with unknown camera constraints.  In International Workshop on Robot Vision, pages 16–28. Springer, 2008.
  * [43] C. Villani.  Topics in optimal transportation.  Number 58. American Mathematical Soc., 2003.
  * [44] C. Villani.  Optimal Transport: Old and New.  Springer, 2008.
  * [45] T. Vogt, R. Haase, D. Bednarski, and J. Lellmann.  On the connection between dynamical optimal transport and functional lifting.  arXiv:2007.02587, 2020.
  * [46] T. Vogt, E. Strekalovskiy, D. Cremers, and J. Lellmann.  Lifting methods for manifold-valued variational problems.  In Handbook of Variational Methods for Nonlinear Geometric Data, pages 95–119. Springer, 2020.
  * [47] M. J. Wainwright, M. I. Jordan, et al.  Graphical models, exponential families, and variational inference.  Foundations and Trends® in Machine Learning, 1(1–2):1–305, 2008.
  * [48] Y. Wald and A. Globerson.  Tightness results for local consistency relaxations in continuous MRFs.  In Conference on Uncertainty in Artificial Intelligence (UAI), 2014\.
  * [49] S. Wang, A. Schwing, and R. Urtasun.  Efficient inference of continuous Markov random fields with polynomial potentials.  In Advances in Neural Information Processing Systems (NeurIPS), 2014\.
  * [50] A. Weinmann, L. Demaret, and M. Storath.  Total variation regularization for manifold-valued data.  SIAM J. Imaging Sci., 7(4):2226–2257, 2014.
  * [51] A. Weinmann, L. Demaret, and M. Storath.  Mumford–Shah and Potts regularization for manifold-valued data.  Journal of Mathematical Imaging and Vision, 55(3):428–445, 2016\.
  * [52] T. Werner.  A linear programming approach to max-sum problem: A review.  IEEE Trans. Pattern Anal. Mach. Intell. (PAMI), 29(7):1165–1179, 2007.
  * [53] K. Yamaguchi, T. Hazan, D. McAllester, and R. Urtasun.  Continuous Markov random fields for robust stereo estimation.  In European Conference on Computer Vision (ECCV), 2012.
  * [54] C. Zach.  Dual decomposition for joint discrete-continuous optimization.  In International Conference on Artificial Intelligence and Statistics (AISTATS), 2013.
  * [55] C. Zach and P. Kohli.  A convex discrete-continuous approach for Markov random fields.  In European Conference on Computer Vision (ECCV), 2012.