Zhenzhang Ye
email: zhenzhang.ye@tum.de

Bjoern Haefner
email: bjoern.haefner@tum.de

Yvain Quéau
email: yvain.queau@ensicaen.fr

Thomas Möllenhoff
email: thomas.moellenhoff@riken.jp

Daniel Cremers
email: cremers@tum.de

A Cutting-Plane Method for Sublabel-Accurate Relaxation of Problems with Product Label Spaces

Keywords: Variational Methods, Manifold-Valued Problems, Convex Relaxation, Global Optimization

 which combined two recent approaches to make lifting more scalable: product-space relaxation and sublabel-accurate discretization. Furthermore, it is shown that a simple cutting-plane method can be used to solve the resulting semi-infinite optimization problem. This journal version extends the previous conference work with additional experiments, a more detailed outline of the complete algorithm and a user-friendly introduction to functional lifting methods.

Introduction

In this paper, we present a convex optimization framework for total-variation regularized problems of the following form: (1)

The set Γ = {(γ 1 , . . . , γ k ) ∈ R N : γ i ∈ Γ i , i = 1 . . . k} is defined by k individual submanifolds

Γ i ⊂ R Ni with N = N 1 + . . . + N k .
The individual Γ i are required to be bounded subsets of R Ni .

Since the focus of this paper are imaging applications we assume Ω ⊂ R 2 to be a rectangular domain but the approach is easily generalized to higher dimensional or non-rectangular domains.

We make no special assumptions on the cost c : Ω×Γ → R ≥0 in (1) and allow it to be a general nonnegative nonconvex function. This turns (1) into an overall nonconvex optimization problem, which can be challenging to solve using standard gradientbased methods. Moreover, we do not assume that we are able to compute gradients, projections or proximal operators of the cost function c(x, u(x)).

Our approach only requires function evaluations. This allows us to consider degenerate costs that are out of reach for gradient-based approaches.

The regularizer in (1) is a separable total variation regularization TV(u i ) on the individual components u i : Ω → R Ni weighted by a parameter λ i > 0. The total variation (TV) [START_REF] Pock | Global solutions of variational models with convex regularization[END_REF][START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] encourages a spatially smooth but edge-preserving solution. It can be defined as the following equation:

TV(u i ) := sup p:Ω→R N i ×2 ∥p(x)∥ * ≤1 Ω ⟨Div x p(x), u i (x)⟩ dx (2) = Ω ∥∇u i (x)∥ dx,
where by ∇u i (x) ∈ R Ni×2 we denote the Jacobian matrix and by ∥ • ∥ * the dual norm of ∥ • ∥. The last equality in (2) holds for sufficiently smooth u i .

The convex relaxation approach we use in this paper works for general convex and nonconvex regularizers which depend on the Jacobian ∇u i , see [START_REF] Möllenhoff | Lifting vectorial variational problems: a natural formulation based on geometric measure theory and discrete exterior calculus[END_REF]; [START_REF] Pock | An algorithm for minimizing the piecewise smooth Mumford-Shah functional[END_REF][START_REF] Pock | Global solutions of variational models with convex regularization[END_REF]; [START_REF] Vogt | Handbook of variational methods for nonlinear geometric data[END_REF]. However, the main focus of this paper is an efficient implementation of the data cost c, and therefore we consider only the separable total variation (2).

Motivation and Applications.

To motivate Problem (1), let us consider some practical applications in low-level vision and imaging. One example is the variational estimation of optical flow between two RGB images I 1 , I 2 : Ω → R 3 , see [START_REF] Horn | Determining optical flow[END_REF]. In that case, Γ 1 = Γ 2 = [a, b] ⊂ R models the displacement between the two images and the cost function is given by a photometric error.

Often, Γ i is a curved manifold, see, e.g., the applications presented by [START_REF] Lellmann | Total variation regularization for functions with values in a manifold[END_REF]; [START_REF] Weinmann | Total variation regularization for manifoldvalued data[END_REF]. Examples include Γ i = S 2 for normal field processing (Lellmann, Strekalovskiy, et al., 2013), SO(3) for motion estimation [START_REF] Görlitz | Piecewise rigid scene flow with implicit motion segmentation[END_REF] or the circle S 1 for processing of cyclic data (Cremers & Strekalovskiy, 2013;[START_REF] Steinke | Nonparametric regression between general Riemannian manifolds[END_REF].

Many real-wolrd applications requires to estimate multiple quantities in a joint fashion. This naturally leads to the formulation of product space which is considered in (1), where Γ = Γ 1 × . . . × Γ k .

In this case, each Γ i models one quantity of interest that one aims to estimate.

A prominent approach to address joint optimization problems of this form are alternating procedures such as expectation maximization [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], block-coordinate descent and alternating direction-type methods [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. There, the idea is to estimates a single quantity while holding the other ones fixed. However, these approaches usually depend on a good initialization and are easy to get stuck in a very poor local optima. Therefore, the goal of this paper is to instead consider a convex relaxation of Problem (1). The relaxed problem can then be solved to global optimality with standard proximal methods such as the primal dual algorithm [START_REF] Pock | Diagonal preconditioning for first order primal-dual algorithms in convex optimization[END_REF]. These methods are usually well parallelizable. Thus, they can be efficiently implemented on GPUs, allowing to solve large-scale problems in a reasonable time.

Contributions.

The main difficulty with convex approaches to (1) is the large memory requirements which are inherent to a lifted problem formulation which renders the problem convex. In order to improve the memory-efficiency of relaxations, two disparate ideas have been considered in previous work: sublabel-accurate liftings [START_REF] Möllenhoff | Sublabel-accurate relaxation of nonconvex energies[END_REF] and product-space relaxations [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF]. In this paper, we combine both approaches and present a sublabel-accurate implementation of [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF]. Unlike previous liftings [START_REF] Möllenhoff | Sublabel-accurate relaxation of nonconvex energies[END_REF][START_REF] Möllenhoff | Sublabelaccurate discretization of nonconvex freediscontinuity problems[END_REF][START_REF] Vogt | Handbook of variational methods for nonlinear geometric data[END_REF], our approach does not require epigraphical projections and can therefore be applied in a black-box fashion, requiring only evaluations of the cost c.

Our main contribution is a simple way to implement the resulting semi-infinite optimization problem with a cutting-plane method. Moreover, we show that using this method, we can achieve a lower energy than the product-space lifting [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF] on optical flow estimation and manifold-valued denoising problems.

This journal paper is an extended version of the conference paper [START_REF] Ye | Sublabel-accurate multilabeling meets product label spaces[END_REF]. In particular, we offer the following contributions over the conference version:

• We have added additional background and explanations on the basics of functional lifting to Section 2 and Section 3. • In Section 4 we added the detailed update equations for the primal-dual algorithm. • We added additional figures and explanations to Section 5 to illustrate and provide an intuition of our algorithm on a simple example.

• We added additional experiments on optical flow and manifold-valued image denoising to Section 5 and evaluated our method on a larger set of images.

Overview of the paper.

After this introduction, we provide in Section 2 an introduction to functional lifting methods for Problem (1), review existing works and explain our contributions relatively to them. We summarize the convex relaxation for (1) and its discretization in Section 3. The proposed cutting-plane method and sampling strategy which we use to implement the discretized relaxation are presented in Section 4. In Section 5 we evaluate our method on a toy problem and several real-world imaging applications. Our conclusions are eventually drawn in Section 6.

An introduction to functional lifting

Let us first consider a simplified version of Problem (1) where Ω consists only of a single point, i.e., the nonconvex minimization of one data term:

min γ∈Γ c(γ 1 , . . . , γ k ). (3)
A well-known approach to the global optimization of (3) is a lifting or stochastic relaxation procedure, which has been considered in diverse fields such as polynomial optimization [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF], continuous Markov random fields [START_REF] Bauermeister | Lifting the convex conjugate in lagrangian relaxations: A tractable approach for continuous Markov random fields[END_REF][START_REF] Fix | Duality and the continuous graphical model[END_REF][START_REF] Peng | Convex max-product algorithms for continuous MRFs with applications to protein folding[END_REF], variational methods [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF], and blackbox optimization [START_REF] Carlier | On a class of multidimensional optimal transportation problems[END_REF][START_REF] Ollivier | Information-geometric optimization algorithms: A unifying picture via invariance principles[END_REF][START_REF] Schaul | Studies in continuous black-box optimization[END_REF]. The idea is to relax the search space in (3) from γ ∈ Γ to probability distributionsu ∈ P(Γ) and solve

1 min u∈P(Γ) Γ c(γ 1 , . . . , γ k ) du(γ 1 , . . . , γ k ). (4)
Due to linearity of the integral wrt. u and convexity of the relaxed search space, this is a convex problem for arbitrary cost c. Moreover, the minimizers of (4) concentrate at the optima of c and can hence be identified with solutions to (3). However, if Γ is a continuum, this problem is infinite-dimensional and therefore challenging.

We illustrate the conceptual difference between the formulation (3) and ( 4) on a one-dimensional example in Figure 1.

Discrete/traditional multilabeling.

In the context of Markov random fields [START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF][START_REF] Kappes | A comparative study of modern inference techniques for discrete energy minimization problems[END_REF] and multilabel optimization (Chambolle et al., 2012;[START_REF] Lellmann | Convex multi-class image labeling by simplex-constrained total variation[END_REF][START_REF] Lellmann | Continuous multiclass labeling approaches and algorithms[END_REF][START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF] one typically discretizes Γ into a finite set of points (called the labels) Γ = {v 1 , . . . , v ℓ }. This turns (4) into a finite-dimensional linear program min u∈∆ ℓ ⟨c ′ , u⟩ where c ′ ∈ R ℓ ≥0 denotes the label cost and ∆ ℓ ⊂ R ℓ is the (ℓ -1)-dimensional unit simplex. If we evaluate the cost at the labels, this program upper bounds the continuous problem (3), since instead of all possible solutions, one considers a restricted subset determined by the labels. Since the solution will be attained at one of the labels, typically a fine meshing is needed. Similar to black-box and zero-order optimization methods, this strategy suffers from the curse of dimensionality. When each Γ i is discretized into ℓ labels, the overall number is ℓ k which quickly becomes intractable since many labels are required for a smooth solution. This limits the method to be applied on more practical problems. Additionally, for pairwise or regularizing terms, often a large number of dual constraints has to be implemented. In that context, the work from Lellmann, Lellmann, et al. (2013) considers a constraint pruning strategy as an offline-preprocessing.

Fig. 1 Traditional optimization vs. optimization over a space of probability distributions. In the top row, going from left to right, we illustrate a traditional gradient-based local optimization method on a nonconvex problem of the form (3). Given a bad initialization, the algorithm might get stuck in a poor local optimum. The bottom row illustrates an optimization procedure on the relaxed problem (4). Due to convexity of the objective and search space of probability measures, the solution concentrates at a Dirac distribution centered around a global optimum.

Sublabel-accurate multilabeling.

The discrete-continuous MRF [START_REF] Fix | Duality and the continuous graphical model[END_REF][START_REF] Zach | Dual decomposition for joint discrete-continuous optimization[END_REF][START_REF] Zach | A convex discretecontinuous approach for Markov random fields[END_REF] and lifting methods [START_REF] Laude | Sublabel-accurate convex relaxation of vectorial multilabel energies[END_REF]Möllenhoff et al., 2016;[START_REF] Möllenhoff | Sublabelaccurate discretization of nonconvex freediscontinuity problems[END_REF] attempt to find a more label-efficient convex formulation. These approaches can be understood through duality [START_REF] Fix | Duality and the continuous graphical model[END_REF][START_REF] Möllenhoff | Sublabelaccurate discretization of nonconvex freediscontinuity problems[END_REF]. Applied to (3), the idea is to replace the cost c : Γ → R with a dual variable q : Γ → R: min u∈P(Γ) sup q:Γ→R Γ q(γ 1 , . . . , γ k ) du(γ 1 , . . . , γ k ), s.t. q(γ) ≤ c(γ) for all γ ∈ Γ.

(5)

The inner supremum in the formulation (5) maximizes the lower-bound q. Additionally, if the dual variable is sufficiently expressive, this problem is actually equivalent to (4). Approximating q, for example with piecewise linear functions on Γ, one arrives at a lower-bound to the nonconvex problem (3). It has been observed in a recent series of works [START_REF] Laude | Sublabel-accurate convex relaxation of vectorial multilabel energies[END_REF][START_REF] Möllenhoff | Sublabel-accurate relaxation of nonconvex energies[END_REF][START_REF] Möllenhoff | Sublabelaccurate discretization of nonconvex freediscontinuity problems[END_REF][START_REF] Vogt | Handbook of variational methods for nonlinear geometric data[END_REF][START_REF] Zach | A convex discretecontinuous approach for Markov random fields[END_REF] that piecewise linear dual variables can lead to smooth solutions even when q (and therefore also u) is defined on a rather coarse mesh. As remarked by [START_REF] Fix | Duality and the continuous graphical model[END_REF]; [START_REF] Laude | Sublabel-accurate convex relaxation of vectorial multilabel energies[END_REF]; [START_REF] Möllenhoff | Sublabel-accurate relaxation of nonconvex energies[END_REF], for an affine dual variable this strategy corresponds to minimizing the convex envelope of the cost, min γ∈Γ c * * (γ), where c * * denotes the Fenchel biconjugate of c.

The implementation of the constraints in (5) can be challenging even in the case of piecewiselinear q. This is partly due to the fact that Problem ( 5) is a semi-infinite optimization problem [START_REF] Blankenship | Infinitely constrained optimization problems[END_REF], i.e., an optimization problem with infinitely many constraints. The works [START_REF] Möllenhoff | Sublabel-accurate relaxation of nonconvex energies[END_REF][START_REF] Zach | A convex discretecontinuous approach for Markov random fields[END_REF] implement the constraints via projections onto the epigraph of the (restricted) conjugate function of the cost within a proximal optimization framework. Such projections are only available in closed form for some choices of c and expensive to compute if the dimension is larger than one [START_REF] Laude | Sublabel-accurate convex relaxation of vectorial multilabel energies[END_REF]. This limits the applicability in a "plug-and-play" fashion.

Product-space liftings.

The product-space lifting approach [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF] attempts to overcome the aforementioned exponential memory requirements of labeling methods in an orthogonal way to the sublabel-based methods. The main idea is to exploit the productspace structure in (1) and optimize over k marginal distributions of the probability measure u ∈ P(Γ), which we denote by u i ∈ P(Γ i ). Applying [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF] to the single data term (3) one arrives at the following relaxation:

min {ui∈P(Γi)} sup {qi:Γi→R} k i=1 Γi q i (γ i ) du i (γ i ) s.t. k i=1 q i (γ i ) ≤ c(γ) for all γ ∈ Γ. (6)
Since one only has to discretize the individual Γ i this substantially reduces the memory require-

ments from O(ℓ N ) to O( k i=1 ℓ Ni ).
While at first glance it seems that the curse of dimensionality is lifted, the difficulty is moved to the dual, where we still have a large (or even infinite) number of constraints. A global implementation of the constraints with Lagrange multipliers as proposed in [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF] again leads to the same exponential dependency on the dimension.

As a side note, readers familiar with optimal transport may notice that the supremum in ( 6) is a multi-marginal transportation problem [START_REF] Carlier | On a class of multidimensional optimal transportation problems[END_REF][START_REF] Villani | Optimal transport: Old and new[END_REF] with transportation cost c. This view is mentioned by [START_REF] Bach | Submodular functions: from discrete to continuous domains[END_REF] where relaxations of form ( 6) are analyzed under submodularity assumptions.

In summary, the sublabel-accurate lifting methods, discrete-continuous MRFs [START_REF] Möllenhoff | Sublabel-accurate relaxation of nonconvex energies[END_REF][START_REF] Zach | A convex discretecontinuous approach for Markov random fields[END_REF] and product-space liftings [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF] all share a common difficulty: implementation of an exponential or even infinite number of constraints on the dual variables.

Summary of our contribution.

Our main contribution is a simple way to implement the dual constraints in an online fashion with a random sampling strategy which we present in Section 4. This allows a black-box implementation, which only requires an evaluation of the cost c and no epigraphical projection operations as in [START_REF] Möllenhoff | Sublabel-accurate relaxation of nonconvex energies[END_REF]; [START_REF] Zach | A convex discretecontinuous approach for Markov random fields[END_REF]. Moreover, the sampling approach allows us to propose and implement a sublabel-accurate variant of the product-space relaxation [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF] which we describe in the following section.

Product-space relaxation

Our starting point is the convex relaxation of (1) presented in [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF]; [START_REF] Strekalovskiy | Convex relaxation of vectorial problems with coupled regularization[END_REF]. In these works, Γ i ⊂ R is chosen to be an interval. We first denote the Lagrangian as:

L({u i }, {q i }, {p i }) = k i=1 Ω Γi q i (x, γ i ) (7) -Div x p i (x, γ i ) du x i (γ i ) dx. (8)
Following [START_REF] Vogt | Handbook of variational methods for nonlinear geometric data[END_REF] we consider a generalization to manifolds Γ i ⊂ R Ni which leads us to the following relaxation:

min {ui:Ω→P(Γi)} sup {qi:Ω×Γi→R} {pi:Ω×Γi→R 2 } L({u i }, {q i }, {p i }), s.t. ∥P Tγ i ∇ γi p i (x, γ i )∥ * ≤ λ i , ∀i, x, γ i , (9)
k i=1 q i (x, γ i ) ≤ c(x, γ), ∀x, γ (10)
This cost function appears similar to (6) explained in the previous section, but there are two differences. First, we now have marginal distributions u i (x) for every x ∈ Ω since we do not consider only a single data term anymore. The notation du x i in (8) denotes the integration against the probability measure u i (x) ∈ P(Γ i ). The variables q i play the same role as in ( 6) and lower-bound the cost under constraint (10). The second difference is the introduction of additional dual variables p i and the term -Div x p i in ( 8). Together with the constraint ( 9), this can be shown to implement the total variation regularization as in [START_REF] Lellmann | Total variation regularization for functions with values in a manifold[END_REF]; [START_REF] Vogt | Handbook of variational methods for nonlinear geometric data[END_REF]. Following [START_REF] Vogt | Handbook of variational methods for nonlinear geometric data[END_REF], the derivative

∇ γi p i (x, γ i ) in (9) denotes the (N i ×2)-dimensional
Jacobian considered in the Euclidean sense and P Tγ i the projection onto the tangent space of Γ i at the point γ i .

To get an intuition on the total variation regularization, we use a concrete example to illustrate how (8) and (9) implement it. Consider the case when the labeling variable u x i = δ u(x) is given as a Dirac measure at every point x. As a concrete example, we consider k = 1 for simplicity. The term in (8) then simplifies to

Ω -Div x p(x, u(x)) dx (11) = Ω ⟨∇ γ p(x, u(x)), ∇u(x)⟩ dx, (12)
which follows by applying the chain-rule and the fact that p has compact support. Finally, taking a point-wise supremum over p inside the integral in ( 12) under the dual-norm constraint (9) gives us the total variation of u: Ω ∥∇u(x)∥ dx, which is the same as defined in (2).

Finite-element discretization

We approximate the infinite-dimensional problem (8) by restricting u i , p i and q i to be piecewise functions on a discrete meshing of Ω × Γ i . The considered discretization is a standard finite-element approach and largely follows the work from [START_REF] Vogt | Handbook of variational methods for nonlinear geometric data[END_REF]. Unlike the forward-differences considered in [START_REF] Vogt | Handbook of variational methods for nonlinear geometric data[END_REF] we use lowest-order Raviart-Thomas elements (see, e.g., [START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF], Section 5) in Ω, which are specifically tailored towards the considered total variation regularization.

Discrete mesh.

We approximate each d i -dimensional manifold Γ i ⊂ R Ni with a simplicial manifold Γ h i , given by the union of a collection of d i -dimensional simplices T i . We denote the number of vertices ("labels") in the triangulation of Γ i as ℓ i . The set of labels is denoted by L i = {v i,1 , . . . , v i,ℓi }. As assumed, Ω ⊂ R 2 is a rectangle which we split into a set of faces F of edge-length h x with edge set E. The number of faces and edges are denoted by F = |F|, E = |E|. Data term and the u i , q i variables.

We assume the cost c : Ω × Γ → R ≥0 is constant in x ∈ Ω on each face and denote its value as c(x(f ), γ) for f ∈ F, where x(f ) ∈ Ω denotes the midpoint of the face f . Similarly, we also assume the variables u i and q i to be constant in x ∈ Ω on each face but continuous piecewise linear functions in γ i . They are represented by coefficient functions u h i , q h i ∈ R F •ℓi , i.e., we specify the values on the labels and linearly interpolate inbetween. This is done by the interpolation operator W i,f,γi : R F •ℓi → R which given an index 1 ≤ i ≤ k, face f , and (continuous) label position γ i ∈ Γ i computes the function value based on barycentric coordinates:

W i,f,γi u h i = u i (x(f ), γ i ). Note that after discretization, u i is only defined on Γ h i but we can uniquely associate to each γ i ∈ Γ h i a point on Γ i .
Divergence and p i variables.

Our variable p i is represented by coefficients p h i ∈ R E•ℓi which live on the edges in Ω and the labels in Γ i . The vector p i (x, γ i ) ∈ R 2 is obtained by linearly interpolating the coefficients on the vertical and horizontal edges of the face and using the interpolated coefficients to evaluate the piecewise-linear function on Γ h i . Under this approximation, the discrete divergence Div h

x : R E

•ℓi → R F •ℓi is given by (Div h x p h i )(f ) = p h i (e r ) + p h i (e t ) -p h i (e l ) -p h i (e b ) /h
x where e r , e t , e l , e b are the right, top, left and bottom edges of f , respectively.

Total variation constraint.

Computing the operator P Tγ i ∇ γi is largely inspired by [START_REF] Vogt | Handbook of variational methods for nonlinear geometric data[END_REF], Section 2.2. It is implemented by a linear map D i,f,α,t : R E•ℓi → R di×2 . Here, f ∈ F and α ∈ [0, 1] 2 correspond to a point x ∈ Ω while t ∈ T i is the simplex containing the point corresponding to γ i ∈ Γ i . First, the operator computes coefficients in R ℓi of two piecewise-linear functions on the manifold by linearly interpolating the values on the edges based on the face index f ∈ F and α ∈ [0, 1] 2 . For each function, the derivative in simplex t ∈ T i on the triangulated manifold is given by the gradient of an affine extension. Projecting the resulting vector onto the d i -dimensional tangent space for both functions leads to a d i × 2-matrix which approximates P Tγ i ∇ γi p i (x, γ i ).

Final discretized problem.

Plugging our discretized u i , q i , p i into (8), we arrive at the following finite-dimensional optimization problem:

min {u h i ∈R F •ℓ i } max {p h i ∈R E•ℓ i }, {q h i ∈R F •ℓ i } h 2 x • k i=1 ⟨u h i , q h i -Div h x p h i ⟩ + f ∈F i{u h i (f ) ∈ ∆ ℓi }, (13) s
.t. ∥D i,f,α,t p h i ∥ * ≤ λ i , ∀i ∈ [k], f ∈ F, α ∈ {0, 1} 2 , t ∈ T i , (14) k i=1 W i,f,γi q h i ≤ c (x(f ), γ) , ∀f ∈ F, γ ∈ Γ, ( 15
)
where i{•} is the indicator function. In our applications, we found that it is sufficient to enforce the constraint (14) at the corners of each face which corresponds to choosing α ∈ {0, 1} 2 . Apart from the infinitely many constraints in (15), this is a finite-dimensional convex-concave saddle-point problem and can be tackled by many numerical optimization algorithms.

Solution recovery

Before presenting in the next section how we propose to implement the constraints (15), we briefly discuss how a primal solution {u h i } of the above problem is turned into an approximate solution to (1). To that end, we follow [START_REF] Lellmann | Total variation regularization for functions with values in a manifold[END_REF]; [START_REF] Vogt | Handbook of variational methods for nonlinear geometric data[END_REF] and compute the Riemannian center of mass via an iteration τ = 1, . . . , T :

V τ j = log u τ i (v i,j ), (16)
v τ = ℓi j=1 u h i (f, j)V τ j , (17)
u τ +1 i = exp u τ i (v τ ). (18)
Here, u 0 i ∈ Γ i is initialized by the label with the highest probability according to u h i (f, •). log u τ i and exp u τ i denote the logarithmic and exponential mapping between Γ h i and its tangent space at u τ i ∈ Γ i , which are both available in closed-form for the manifolds we consider here. In our case T = 20 was enough to reach convergence. For flat manifolds, T = 1 is enough, as both mappings boil down to the identity and ( 18) computes a weighted Euclidean mean.

In general, there is no theory which shows that u T (x) = (u T 1 (x), . . . , u T k (x)) from ( 18) is a global minimizer of (1). Tightness of the relaxation in the special case k = 1 and Γ ⊂ R is shown in [START_REF] Pock | Global solutions of variational models with convex regularization[END_REF]. For higher dimensional Γ, the tightness of related relaxations is ongoing research; see [START_REF] Ghoussoub | Hidden convexity in a problem of nonlinear elasticity[END_REF] for results on the Dirichlet energy. By computing a-posteriori optimality gaps, solutions of (8) were shown to be typically near the global optimum of Problem (1); see, e.g., [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF].

Implementation of the constraints

Though the optimization variables in (13) are finitedimensional, the energy is still difficult to optimize because of the infinitely many constraints in (15). Before we present our approach, let us first describe what we refer to as the baseline method for the remainder of this paper. As the baseline approach, we consider the direct solution of (13) where we implemented the constraints only at the label/discretization points L 1 × . . . × L k via Lagrange multipliers (this strategy is also employed by the global variant of the productspace approach [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF]). This baseline actually corresponds to a single outer iteration N it of the proposed Algorithm 1, with a large number M it of inner iterations.

We aim for a framework that allows for solving a better approximation of (15) than the baseline while being of similar memory complexity. To this end, Algorithm 1 alternates the following two steps:

1) Sampling. Based on the current solution we prune previously considered but feasible constraints and sample a new subset of the infinitely many constraints in (15). From all the current sampled constraints, we consider the most violated constraints for each face, add one sample at the current solution and discard the rest.

2) Solving the subsampled problem. Considering the current finite subset of constraints, we solve Problem (13) using a primal-dual method (Chambolle & Pock, 2011) with diagonal preconditioning [START_REF] Pock | Diagonal preconditioning for first order primal-dual algorithms in convex optimization[END_REF]. Both constraints ( 14) and ( 15) are implemented using Lagrange multipliers.

These two phases are performed alternatingly, with the aim to eventually approach the solution of the continuous problem (13). In practice, a fixed number of outer iterations N it is set. While we do not prove convergence of the overall algorithm, convergence results for related procedures exist; see, e.g., [START_REF] Blankenship | Infinitely constrained optimization problems[END_REF], Theorem 2.4.

The detailed algorithm is explained in Algorithm 1. The cost matrix C is constructed by evaluating c(x(f ), γ) at proposed samples S f . We denote ξ and ν as the Lagrange multipliers. The Lagrange multiplier ξ is initialized by a warm-start strategy, i.e. ξ it keeps same if we have the same proposed sample from previous outer iteraion. The prox of a function g with step size τ is defined as:

prox τ g (x) = arg min y 1 2τ ∥x -y∥ + g(y) (19)
Our constraint sampling strategy is detailed in Algorithm 2. For each face in F, it generates a finite set of "sublabels" S f ⊂ Γ at which we implement the constraints (15). Next, we provide the motivation behind each line in the algorithm.

Algorithm 1 Proposed algorithm for problem (1).

Input:

c : Ω × Γ → R, λ i > 0, N it > 0, M it > 0, n > 0, δ > 0, r > 0. 1: u h,0 i = 1/ℓ i , q h,0 i = 0, p h,0 i = 0. 2: S 0 f = L 1 × . . . × L k . 3: for it = 0 to N it do 4:
Construct interpolation matrix W it i and cost matrix C it based on S it f from Algorithm 2.

5:

Initialize ξ it i with the warm-start strategy.

6:

Construct the diagonal preconditioners T it u , T it ν , T it ξ , Σ it p and Σ it q by Pock and Chambolle (2011).

7:

for j = 0 to M it do 8:

p h,j+1 i = p h,j i + Σ it p (-h 2 x • (Div h x ) T u h,j i + D T i,α,t ν j i ) 9: q h,j+1 i = q h,j i + Σ it q (h 2 x • u h,j i -(W it i ) T ξ j i ) 10: ph,j i = 2p h,j+1 i -p h,j i 11: qh,j i = 2q h,j+1 i -q h,j i 12: u h,j+1 i = prox •∈∆ ℓ i (u h,j i -T it u (q h,j i -Div h x ph,j i ))
13:

ν j+1 i = prox λT it ν ∥•∥2 (ν j i -T it ν D i,v,t ph,j i )
14:

ξ j+1 i = prox •≥0 (ξ j i -T it ξ (C it - k i W it i qh,j i ))
15:

end for 16:

Get sampled S it+1 f for each face by Algorithm 2. 17: end for Algorithm 2 Sampling strategy at face f ∈ F.

Input:

Solution u = (u 1 , . . . , u k ) at face f , sublabel-set S f , n, δ, r. /* global exploration */ 1: S ′ f ← uniformSample(Γ, n) /* local exploration around solution */ 2: S ′ f ← S ′ f ∪ localPerturb(u, δ, n) /* remove feasible constraints */ 3: S f ← γ ∈ S f : k i=1 q i (f, γ) > c(f, γ) /* add the most violated r samples */ 4: S f ← top-k(S ′ f , r) ∪ S f /* have one sample at current solution */ 5: S f ← S f ∪ {u}. 6: Return S f

Random uniform sampling (Line 1).

To have a global view of the cost function, we consider a uniform sampling on the label space Γ. The parameter n > 0 determines the number of the samples for each face.

Local perturbation around the mean (Line 2).

Besides the global information, we apply local perturbation around the current solution u. In case the current solution is close to the optimal one, this strategy allows us to refine it with these samples. The parameter δ > 0 determines the size of the local neighbourhood. In our experiments, we always used a Gaussian perturbation with δ = 0.1.

Pruning strategy (Lines 3-4).

Most samples from previous iterations are discarded because the corresponding constraints are already satisfied. We prune all current feasible constraints as in [START_REF] Blankenship | Infinitely constrained optimization problems[END_REF]. Similarly, the two random sampling strategies (Lines 1 and 2) might return some samples for which the constraints are already fulfilled. Therefore, we only consider the samples with violated constraints and pick the r most violated from them. This pruning strategy is essential for a memory efficient implementation as shown later.

Sampling at u (Line 5).

Finally, we add one sample which is exactly at the current solution u ∈ Γ to have at least one guaranteed sample per face. In the next section, we illustrate the behavior of Algorithm 1 on a toy problem, and evaluate its performance on realworld imaging problems. The dual variable q satisfies the constraints on the samples. However, q is not globally optimal as it violates the constraint on the optimal solution ("green > blue").

Numerical validation

Our approach and the baseline are implemented in PyTorch. Code for reproducing the following experiments can be found here: https://github.com/ zhenzhangye/sublabel_meets_product_space. Note that that one of the runtime bottlenecks of our sampling strategy is creating the samples and picking the most violated r as shown in Algorithm 2. Additionally, the sparse matrix operations and the PDHG updates can be more efficiently implemented in CUDA, as all PDHG updates can be executed in a single CUDA kernel, compared to PyTorch's multiple kernel calls. Therefore, a specialized implementation as in [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF] will allow the method to scale by factor 10 -100× in favor of runtime.

Illustration of the Algorithm

First of all, we consider a simplistic minimization problem on a single nonconvex data term:

c(u) = min      -4u + 2.4, u ∈ [0.1, 0.35), 4u -0.4, u ∈ [0.35, 0.6], 2, otherwise, (20)
with 5 labels to illustrate the behavior of both, the baseline algorithm and our sampling strategies. Figure 2 depicts the baseline's behavior. While it only evaluates the energy on the labels, five samples are considered as illustrated by the red dots in Figure 2 (a). Figure 2 (b) shows the dual variable q h i after it iterations. Since the algorithm is maximizing q h i , the green and red dots should overlay (e.g. the second label) when it converges. Despite the convergence, the resulting q h violates the constraints significantly close to the optimal solution ("green > blue"). Therefore, to attain the global optimal solution, the baseline approach needs more labels which requires more memory.

The motivation of random uniform sampling (Line 1, Algorithm 2) and local perturbation around the mean (Line 2, Algorithm 2) in our sampling strategy is intuitively clear. However, as demonstrated later in the experiment of a truncated quadratic energy, sampling at u (Line 5, Algorithm 2) is critical for the stability of our method. A comparison in Figure 3 helps to show the necessity of this strategy. We ran two experiments with the identical settings except for the sampling at u. After a given number of iterations, the dual variable q h is approximately optimal, as indicated in Figure 3 (a). Our pruning strategy (Line 3, Algorithm 2) removes all of the proposed samples since all of them satisfy the constraints. As a result, the subproblem becomes unconstrainted on that dual variable q h and its update has no significance, Figure 3 (b). To solve this problem, we propose to always at least have one sample at u even when q h is nearly optimal, cf. Figure 3 (c). As illustrated in Figure 3 (d), this can avoid the degeneration of q h as it is still constrained.

Finally, the complete sampling strategy is illustrated in Figure 4. As shown in Figure 4 (a), the primal-dual method can obtain the optimal q h for the sampled subproblem. Our sampling strategy can provide necessary samples and prune the feasible ones, cf. Figure 4 (b). These few but meaningful samples lead the q h to achieve global optimality, cf. Figure 4 (c).

Truncated Quadratic Energy

In this section, we study the numerical effect of each line in Algorithm 2. We evaluate our method on the truncated quadratic energy c(x, u(x)) = min{(u(x) -f (x)) 2 , ν}. where f : Ω → R is the input data as show in Figure 5. For this specific experiment, we generate a 64 × 64 gray image degraded with Gaussian noise of standard deviation σ = 0.05 and 5% salt-and-pepper noise. The parameters are chosen as ν = 0.025, λ = 0.25, N it = 10, M it = 200, n = 10 and r = 1. To reduce the effect of randomness, we run each algorithm 20 times and report mean and standard deviation of the final energy for different number of labels in Table 1. We want to emphasize that more labels have benefits for both baseline and our algorithm.

Nevertheless, the proposed approach can reach lower energies with the same number of labels and similar memory requirements.

As can be seen in this table, adding uniform sampling and picking the most violated constraint per face (Lines 1 and 4 of Algorithm 2) already decreases the final energy significantly. We also consider local exploration around the current solution (Line 2), which helps to find better energies at the expense of higher memory requirements.

The pruning strategy (Line 3) circumvents this memory issue, however the energy deteriorates dramatically because some faces could end up having no samples after pruning. Therefore, keeping the current solution as a sample (Line 5) per face prevents the energy from degrading. Including all these sampling strategies, the proposed method can achieve the best energy and runtime, at comparable memory usage to the baseline method.

We further illustrate the comparison on the number of iterations and time between the baseline and our proposed method in Figure 6. Due to the replacement on the samples, we have a peak right after each sampling phase. The energy however converges immediately, leading to an overall decreasing trend.

Additionally, we compare our method to the baseline on a more practical dataset CBSD from [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. This dataset contains 68 images and noisy ones with additive white Gaussian noise (5% in this experiment). The number of labels is 7 for both methods. 5K iterations are performed on the baseline method, while we set N it = 10 and M it = 500 to get a fair comparison. The other parameters are chosen as λ = 0.25, n = 50 and r = 1. The results are shown in Figure 7. Our method outpeforms the baseline among all the images regarding both energy and peak signal-to-noise ratio (PSNR).

Manifold-value Denoising

To show the flexibility of our algorithm, we next evaluate it on a manifold-valued denoising problem in HSV color space. The hue component of this space is a circle, i.e.,

Γ 1 = S 1 , Γ 2 , Γ 3 = [0, 1].
The data term of this experiment is still a truncated quadratic energy, where for the hue component the distance is taken on the circle S 1 . The input images [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF][START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] are degraded with the same setting as above.

Both the baseline and our method are implemented with 7 labels. First of all, we evaluate the impact of the most violated r samples. As shown in Figure 8, The maximum difference of energy and memory is only 0.8% and 0.06%, respectively, which can be considered almost constant wrt. r. Therefore, we pick r = 5. To get an equal number of total iterations, 30K iterations are performed on the baseline, while we set N it = 100 outer iterations with M it = 300 inner primal-dual steps for our method. Other parameters are chosen as λ = 0.015 and n = 30. As shown in Figure 9, our method can achieve a lower energy than the baseline. Qualitatively, since our method implements the constraints not only at the labels but also inbetween, there is less bias.

Optical flow

Given two input images I 1 , I 2 , we compute the optical flow u : Ω → R 2 for the label space Γ = [a, b] 2 . We use a simple ℓ 1 -norm for the data term, i.e. c(x, u(x)) = ∥I 2 (x) -I 1 (x + u(x))∥ 1 and set the regularization weight as λ = 0.04. The baseline approach runs for 50K iterations, while we set N it = 50 and M it = 1000 for a fair comparison. Additionally, the parameters are chosen as n = 20 and r = 1 in Algorithm 2.

We consider three methods for this experiment: the baseline, the method from [START_REF] Laude | Sublabel-accurate convex relaxation of vectorial multilabel energies[END_REF] and ours. The method from [START_REF] Laude | Sublabel-accurate convex relaxation of vectorial multilabel energies[END_REF] approximates the dataterm in a piecewise convex manner and requires a specific epigraph projection for the dataterm. Though it can attain lower energy, our approach requires less memory and tackle any cost function in a simple black-box fashion. Table 2 summarizes the quantitative results obtained on the Middleburry dataset (Baker et Fig. 6 Comparison between the baseline and our approach on a 64 × 64 grayscale image shown in Figure 5, degraded with Gaussian and salt-and-pepper noise. Our approach finds lower energies in fewer iterations and less time than the baseline, which implements the constraints only at the label points. Table 2 We compute the optical flow on the Middlebury dataset [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF] using our method and the baseline for a varying amount of labels. Given an equal number of labels/memory, our sampling strategy performs favorably to an implementation of the constraints at the labels. The relative numbers of energy, memory, average endpoint error (aep) and average angular error (aae) are calculated as "mean( Ours Baseline )" across all 8 datasets. The number in the parentheses resemble the standard deviation. The detailed table with all results can be found in the appendix.

2011), while the detailed absolute numbers can be found in the appendix. This table shows how our approach performs relatively to the baseline, i.e. "mean( Ours Baseline )", e.g. for three labels our energy is 50.91% of the baseline energy for all 8 datasets, while using 99.84% of the baseline's memory and having an average end point error (aep) and average angular error (aae) of 91.94% and 82.34% of the baseline error metrics, respectively. To enable qualitative comparison, we visualize in Figure 10 the results on two of the datasets. The remaining qualitative results on the Middlebury data set [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF] are shown in the appendix. Our method outperforms the baseline approach regarding energy under the same number of labels and requires the same amount of memory. Because [START_REF] Laude | Sublabel-accurate convex relaxation of vectorial multilabel energies[END_REF] uses a tighter relxation on the label space, they can achieve lower energy with a smaller number of labels. However, it runs out of memory easily while the proposed method scales better wrt. memory consumption. PSNR Difference (%) Fig. 7 Quantitative results on the CBSD dataset [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. We ordered both results from the best to the worst. The difference is calculated using the formula ours -baseline baseline

. It is clear that our method consistently outperforms the baseline across all images.

Conclusion and Limitations

In this paper we made functional lifting methods more scalable by combining two advances, namely product-space relaxations [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF] and sublabel-accurate discretizations [START_REF] Möllenhoff | Sublabelaccurate discretization of nonconvex freediscontinuity problems[END_REF][START_REF] Vogt | Handbook of variational methods for nonlinear geometric data[END_REF]. This combination is enabled by adapting a cutting-plane method from semi-infinite programming [START_REF] Blankenship | Infinitely constrained optimization problems[END_REF]. This allows an implementation of sublabel-accurate methods without difficult epigraphical projections. Moreover, our approach makes sublabelaccurate functional-lifting methods applicable to any cost function in a simple black-box fashion. In experiments, we demonstrate the effectiveness of the approach over a baseline based on the productspace relaxation [START_REF] Goldluecke | Tight convex relaxations for vectorvalued labeling[END_REF] and provided a proof-of-concept experiment showcasing the method in the manifold-valued setting.

Future work will concentrate on applying and adapting the presented framework to solve large inverse problems in computer vision with multiple data terms, different regularizers and several manifold-valued optimization variables in a joint fashion. However, it is not obvious if our presented cutting plane approach is easily applicable for such large problems or if novel ideas have to be pursued.

A Additional Optical Flow Results

Table 3 shows the energy and memory requirement as well as the average endpoint error (aep) and average angular error (aae) for the baseline, the method from [START_REF] Laude | Sublabel-accurate convex relaxation of vectorial multilabel energies[END_REF] and our approach across different number of labels for all 8 used Middlebury datasets [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF]. Figure 11 and Figure 12 visualize the remaining optical flow results on the Middlebury dataset [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF] using the baseline, the method from Laude et al. ( 2016) and ours for a varying amount of labels. Note that the approach from Laude et al. ( 2016) leverages a tighter discretization on the label space and is implemented on CUDA. Though their approach achieves better energy under fewer labels, ours has a better scability. Table 3 We compute the optical flow on the Middlebury dataset [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF] using the baseline, the method from Laude et al. ( 2016) and ours for a varying amount of labels. We denote the out of memory error as OOM. Given an equal number of labels/memory, our sampling strategy performs favorably to an implementation of the constraints at the labels comparing to the baseline. Additionally, our method obtains a better scalibilty than the one from [START_REF] Laude | Sublabel-accurate convex relaxation of vectorial multilabel energies[END_REF].

Fig. 2

 2 Fig.2Illustration of the baseline algorithm. (a) Five samples (red dots) from the labels are considered (b) The dual variable q satisfies the constraints on the samples. However, q is not globally optimal as it violates the constraint on the optimal solution ("green > blue").

Fig. 3 Fig. 4

 34 Fig.3Benefit of sampling at u (Line 5, Algorithm 2). (a) Because all of the proposed samples (gray dots) fullfil the constraint, they are all pruned (Line 3, Algorithm 2). (b) The updated q deteriorates because the subproblem is unconstrainted on it. (c) At least one sample is taken (namely it , red dot). (d) This sample constraints q and thus prevents a degenerate solution.

Fig. 5

 5 Fig. 5 (a) The original 64 × 64 grayscale image. (b) Degraded with Gaussian noise of standard deviation σ = 0.05 and 5% salt-and-pepper noise.

Fig. 8

 8 Fig.8Comparison between different number of most violated constraints on the energy and memory: (a) the change of the energy with r = 1 as the basis; (b) the change of memory with r = 1 as the basis. It can be observed that both the energy and memory vary very little (0.8% and 0.06%, respectively) regarding different r.

Fig. 10

 10 Fig. 10 Visualization of the optical flow on Grove3 and RubberWhale from the Middleburry dataset (Baker et al., 2011), using the baseline, the method from Laude et al. (2016) and our method for a varying amount of labels. OOM stands for out of memory. More qualitative results can be found in the appendix.

Fig. 11

 11 Fig. 11 Part I: We visualize the optical flow on the Middlebury dataset (Baker et al., 2011) using baseline, the method from Laude et al. (2016) and ours for a varying amount of labels for qualitative inspection. OOM stands for out of memory.

Table 1

 1 Ablation study indicating the effect of individual lines in Algorithm 2. Numbers in parentheses indicate the standard deviation across 20 runs.

		Labels	Baseline	+Line 1&4	+Line 2	+Line 3	+Line 5
	Energy		4589 (±0.00) 2305 (±3.73) 2291 (±3.6) 8585 (±130.4) 2051 (±10.7)
	Time [s]	3	8.98	22.77	23.22	23.22	23.33
	Mem. [Mb]		11.21	13.94	15.53	11.65	12.05
	Energy		2582 (±0.00) 2020 (±2.68) 2012 (±1.3) 7209 (±116.7)	1969 (±3.6)
	Time [s]	7	74.13	16.02	16.61	15.56	18.38
	Mem. [Mb]		28.35	32.96	33.49	28.356	28.68
	Energy		2029 (±0.00) 1935 (±1.14) 1926 (±0.7)	5976 (±75.7)	1901 (±3.7)
	Time [s]	13	183.80	37.65	38.84	38.29	38.22
	Mem. [Mb]		52.85	60.55	60.94	54.35	54.73

P(Γ) is the set of nonnegative Radon measures on Γ with total mass u(Γ) = 1.

Acknowledgements This work was supported by the ERC Advanced Grant SIMULACRON.

Declarations

Conflict of interest The authors declare that they have no conflict of interest. [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF][START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] in HSV color space (Γ 1 = S 1 , Γ 2 = Γ 3 = [0, 1]) using our method and the baseline with 7 labels. The hue component plot demonstrates that our approach is able to handle the manifold setting where the jump from 2π to 0 is permitted. Since our approach implements the constraints adaptively inbetween the labels (gray lines) it reaches a lower energy with less label bias.

Anal., 10 (2), 517-530.

Chambolle, A., Caselles, V., Cremers, D., Novaga, M., [START_REF] Pock | Global solutions of variational models with convex regularization[END_REF]. An introduction to total variation for image analysis. Theoretical Foundations and Numerical Methods for Sparse Recovery, 9 (263-340), 227.

Chambolle, A., Cremers, D., Pock, T. (2012).

A convex approach to minimal partitions.