Composite Optimization by Nonconvex
Majorization-Minimization

arXiv:1802.07072v1 [math.OC] 20 Feb 2018

Jonas Geiping∗

Michael Moeller†

February 21, 2018

Abstract
Many tasks in imaging can be modeled via the minimization of a nonconvex composite function. A popular class of algorithms for solving such problems are majorizationminimization techniques which iteratively approximate the composite nonconvex function by a majorizing function that is easy to minimize. Most techniques, e.g. gradient descent, utilize convex majorizers in order guarantee that the majorizer is easy to
minimize. In our work we consider a natural class of nonconvex majorizers for these
functions, and show that these majorizers are still sufficient for a provably convergent
optimization scheme under quite general assumptions. Numerical results illustrate that
by applying this scheme, one can often obtain superior local optima compared to previous
majorization-minimization methods, when the nonconvex majorizers are solved to global
optimality. Finally, we illustrate the behavior of our algorithm for depth super-resolution
from raw time-of-flight data.

Keywords: Nonconvex Optimization, First-Order Optimization, Majorization-Minimization,
Kurdyka Lojasiewicz inequality, Time-of-Flight Depth Reconstruction
AMS Subject Classification: 90C26, 90C06, 68U10, 32B20, 65K10, 47J06

1

Introduction

Nonlinear optimization methods are a staple for many imaging tasks that can be regarded
as the minimization of some objective function, also called energy. Unfortunately, many
energies arising from the faithful modeling of the data formation process and a state-ofthe-art regularization term are inherently nonconvex, coupled, and high dimensional. Since
determining the global minimizer of such cost function is rarely feasible, one frequently turns
to (gradient-based) methods that only find a (possibly sub-optimal) critical point of the
energy landscape [43].
Interestingly, some high-dimensional nonconvex optimization problems do admit a global
solution within reasonable time. Besides problems for which the solution can be determined analytically, the aforementioned class includes separable problems on a bounded domain, i.e. problems for which the minimization of an energy E with respect to some variable u ∈ Rn decomposes into the minimization of separate low-dimensional energies, e.g.
∗

University of Siegen, Hölderlinstraße 3,
michael.moeller@uni-siegen.de)

57076 Siegen,

1

Germany (jonas.geiping@uni-siegen.de,

(a)

(b)

Figure 1: Nonconvex versus convex majorization. (a) shows an energy of type (2) with a
convex majorizer. (b) shows the same energy, but with a solvable nonconvex majorizer. The
initial point is marked in red, the global minimum of the energy in green. We can see that
the shown nonconvex majorizer can better represent the given function.
P
E(u) = ni=1 Ei (ui ). Even more remarkably, there are several types of non-separable nonconvex optimization problems which can be reformulated as convex problems, e.g. via convex
relaxation techniques [15] or via functional lifting [50], and still yield a globally optimal solution to the original nonconvex problem. Unfortunately, the aforementioned techniques rely
on a special structure of the objective. Even seemingly minor perturbations of the required
structure make it impossible to exploit these techniques, and lead practitioners to consider
local (gradient-based) methods again.
Interestingly, many of such local methods admit an interpretation in the framework of
majorization-minimization techniques: In each iteration, the energy E is approximated by a
simpler function Euk which satisfies
Euk (u) ≥ E(u),
Euk (uk ) = E(uk ),
for uk being the current iterate. By defining the next iterate to be the minimizer of the
approximation Euk ,
uk+1 = arg min Euk (u),
u

one automatically obtains monotonically decreasing objective values.
Common gradient-based methods use simple convex approximation functions Euk , e.g.
quadratic functions,
Euk (u) = E(uk ) + h∇E(uk ), u − uk i +
in the case of gradient descent.
imation functions Euk are only
information about the shape of
In this work we propose a
functions Euk with the idea to

1
||u − uk ||2 ,
2τ

(1)

While this leads to easy-to-solve subproblems, such approxa crude approximation of the original energy and almost all
the original energy landscape is lost.
novel majorization-minimization technique with nonconvex

1. approximate the original energy landscape much more faithfully, and
2

2. still be able to minimize Euk globally by considering functions Euk that are either
separable or can be minimized via relaxation techniques.
As illustrated in a simple two-dimensional example in Figure 1, one can expect a more
faithful approximation of the original energy to yield ’better’ local minima: While the parabolic approximation of Figure 1a yields a nearby local minimum, the separable nonconvex
majorizer in Figure 1b allows to skip several local minima. In this example, the minimizer of
the nonconvex majorizer is in a close vicinity to the global minimizer after just a single step
of the algorithm.
While our motivation comes from the above (somewhat heuristic) idea of finding ’better’
local minima, our convergence analysis does not depend on the subproblems being solved to
global optimality. For the remainder of the paper we consider the minimization of composite
energies of the form
E(u) = G(ρ(u)) + R(u),
(2)
for suitable functions G : Rm → R, ρ : Rn → Rm and R : Rn → R, via the iterative
minimization of
Euk (u) = G(ρ(uk )) + h∇G(ρ(uk ), ρ(u) − ρ(uk )i + R(u) +

1
||ρ(u) − ρ(uk )||2 .
2τ

(3)

The model function Euk is a naturally global, but nonconvex, majorizer of E for suitable τ
as we will see later. A typical example for ’simple’ functions ρ : Rn → Rn and R : Rn → R
is
Pngiven when both functions are separable, i.e. ρ(u) = (ρ1 (u1 ), . . . , ρn (un )) and R(u) =
i=1 ri (ui ). In this case, the nonconvex majorizer (3) is then also separable and can be
solved in each dimension separately.
We continue summarizing some of the related work for nonconvex and composite optimization problems and illustrate how the proposed majorization-minimization technique (3)
differs from the methods that have been considered in the literature so far.

1.1

Related Work

The current field of nonlinear optimization is quite wide. In the following overview of related
work we focus on results, that like our method do not require convexity of the objective
function and we limit us to generalizations of first-order methods. The general framework of
majorization-minimization methods has been reviewed widely in the literature of the recent
decades, see for example, [24, 39, 57, 61].
The first option for tackling the minimization of (2) is to ignore the composite structure
of G ◦ ρ, naturally leading to schemes like the aforementioned gradient descent (GD) (1)
or the closely related forward backward splitting (FBS) [16, 6, 44]. As we will see in more
detail below, the proposed scheme recovers such algorithms in the special case of ρ being
the identity. The convergence1 of a general class of nonconvex first-order descent methods,
including GD and FBS, was shown e.g. in [2]. It is important to note that such a convergence
is nontrivial for arbitrary nonconvex functions and requires, for example, some algebraic
notion of ’tameness’ [27], that is nevertheless usually present in practice.
1

In the context of first order methods, we always consider ’convergence’ as implying that the sequence of
iterates converges to a stationary point of the objective function.

3

The most limiting assumption in these first-order methods is the Lipschitz continuity of
the gradient of F = G ◦ ρ, the first part of the objective function. This class of problems
was recently extended in [9] to L-smooth adaptable function [5], these functions are not
necessarily convex or L-smooth, only a smooth and convex ’kernel’ function h must exist,
so that Lh − F is convex for some L > 0. The previously mentioned methods can be
extended to a descent relative to these kernel functions. Defining the Bregman distance of h
as Dh (u, v) = h(u) − h(v) − h∇h, u − vi, [9]’s majorizer can be written as
1
Euk (u) = h∇F (uk ), u − uk i + R(u) + Dh (u, uk ).
τ

(4)

They show that the sequence of iterates generated by this type of majorizer converges for
appropriate τ and conditions to h, F and R. This is a quite satisfying result as this formulation
incorporates logarithmic or quadric data terms and vector-valued step sizes.
Recent works have also proposed general frameworks for iteratively replacing the original minimization problem with simple approximation functions Euk beyond majorizationminimization. [18] analyzed approximation functions Euk , satisfying |Euk (u)−E(u)| ≤ ω(||u−
uk ||) for a proper growth function ω. A minimization scheme of these approximation functions
exhibits subsequential convergence to critical points, even if the subproblem evaluations are
inexact. These approximation functions need not necessarily be convex, but the distance of
their subsequent evaluations must tend to zero. A slightly different generalization is discussed
in [49], where approximation functions satisfying |Euk (u) − E(u) − Dh (u, uk )| ≤ ω(||u − uk ||)
are examined. Here ω is a growth function and Dh a Bregman distance generated by a ’kernel’ function similarly to the previously discussed (4). Subsequence convergence can again
be shown here, under relatively weak conditions. However the approximation function Euk
is taken to be convex in this work to, among other properties, guarantee the success of a
backtracking scheme and reach an implementable algorithm.
Coming to related work in composite optimization we find that there are two ways to
handle problems of type (2): Either we linearize the outer function G in each approximation,
or the inner function ρ. Linearizing the inner function ρ leads to methods that are reminiscent of classical Levenberg-Marquardt algorithms for nonlinear least-squares problems. The
approximation function can be written as


1
Euk (u) = G ρ(uk ) + Jρ (uk )(u − uk ) + R(u) + ||u − uk ||2 ,
(5)
2τ
where Jρ denotes the Jacobian of ρ. A classical application for this composition are systems
of nonlinear equations. Due to the inner linearization, it is in general not required that G
is smooth. Subsequence convergence follows as a result of [36, 18] or [49]. Further literature
can be found under the terms ’prox-linear’ or ’prox-descent’, e.g. [36, 19]. Linearizing the
outer function leads to algorithms related to iterative re-weighting procedures:
Euk (u) = G(ρ(uk )) + h∇G(ρ(uk )), ρ(u) − ρ(uk ) + R(u) +

1
||u − uk ||2 .
2τ

(6)

Subsequence convergence follows from the general result of [18] under the assumption that
the distance of subsequent iterates tends to zero. Further analysis, related to special cases in
iterative re-weighting can be found in [48] or under more general assumptions, but including
the convexity of Euk in [49]. The connection to iterative reweighting is immediate for concave
4

(a)

(b)

(c)

Figure 2: Visualization of related work. (a) shows the original function of type (2), (b) shows
a gradient descent majorizer (1), (c) shows a forward-backward splitting majorizer (4). The
point uk is equal in each figure and shown in red and the minimizer uk+1 in green.

(a)

(b)

(c)

Figure 3: Visualization of related work. (a) shows a prox-linear type inner linearization
(5), (b) shows an outer linearization (6), (c) shows finally shows the proposed (separable)
majorizer (3). The point uk is equal in each figure and shown in red and the minimizer uk+1
in green.
G, as then τ can be taken arbitrarily large and the proximal term vanishes. This formulation
is closely related to our work and differs from ours in the way we measure the distance to the
previous iterate. We later discuss the implications of this difference.
As a first intuitive visualization, Figure 2 and Figure 3 show these majorization functions
in two dimensions. For a nonconvex function of type (2) in Figure 2a, a gradient descent
majorizer is shown in Figure 2b and a forward-backward splitting in Figure 2c. We clearly
see that for both majorizers their respective minimizers, marked in green, are located in a
close neighborhood to the current iterate, marked in red. Both algorithms will likely converge
to a nearby local minimum of the original energy Figure 2a.
The presented related majorizers for composite optimization are shown in Figure 3. Figure 3a and Figure 3b show both linearization variants, namely (5) and (6). These generally
produce more faithful representations of the original energy (Figure 2a), but both minimizers
are still far away from the global minimum. Finally, Figure 3c shows our majorizer (3). Note

5

that the minimizer of this majorizer can not only be computed efficiently due to its separability, but also allows for a global view of the function and its minimizer almost coincides
with the global minimum although the initial point is quite far from it.
Finally, a recent preprint [8] proposes to solve composite minimization problems with a
different approach, namely a nonlinear splitting variant, reformulating the problem to
min G(v) + R(u)

u,v∈Rn

s.t. ρ(u) = v,

(7)

and introducing an augmented Lagrangian formulation
τ
minn G(v) + R(u) + hw, ρ(u) − vi + ||ρ(u) − v||2 .
u,v∈R
2

(8)

with an additional variable w ∈ Rn that mimics the dual variable of the convex setting. This is
a quite interesting result, as it shows that the complementarity of forward-backward splitting
and augmented Lagrangian methods can be extended into the composite setting. Whereas
our method is a generalization of forward-backward splittings, their work generalizes ADMM
[21]. Critically both ours and their approach rely on the efficient solution of a nonlinear
programming task as intermediate step in the algorithm. For us, this is the nonconvex
majorizer (3), the corresponding problem in [8, Eq. (6.3)] is the minimization of (8) for u:
uk+1 = arg min G(v k+1 ) + R(u) + hwk , ρ(u) − v k+1 i +
u

τk
µ
||ρ(u) − v k+1 ||2 + ||u − uk ||2 . (9)
2
2

Both subproblems are in general equally difficult as they are connected for µ = 0, identifying
v k+1 = ρ(uk ) and wk = ∇G(ρ(uk )).
Although formulated in less generality in terms of the involved functions but in more
generality in terms of the number of nested functions, the update equation of the related
work [20, Eq. (11)] for solving problem (2) can be written as
1
1
uk+1 = arg min R(u) + kρ(u) − ρ(uk )k2 + σh∇G(ρ(uk )), ρ(u) − ρ(uk )i + ku − uk k2 ,
u
2
2τ
for an affine linear ρ. This is similar to the proposed algorithm but also contains the additional
proximity term for u − uk . The analysis we provide in this work could make it interesting to
revisit [20] in the two-layer case.
Solving the subproblems via lifting. While the convergence analysis of our approach
will make rather weak assumptions on the global quality of the solution used in each of
the subproblems (3), we found our method to be particularly effective and successful if the
(nonconvex) subproblems are solved to global optimality. This raises the question what types
of functions allow to determine globally optimal solutions to such subproblems.
A rather simple case occurs if the involved functions are separable or separable into blocks
of few variables. In these situations we can apply exhaustive search and branch-and bound
algorithms to each block separately [29, 23].
More interesting for imaging tasks is the class of functions where the subproblems can be
solved by functional lifting. It was shown in [50, 13] that free discontinuity-type energies, in
particular,
Z
ν(x, u(x), ∇u(x)) dx,

E(u) =
Ω

6

(10)

with ν : R×R×Rd → R being continuous in its second argument, and convex and continuous
in its third argument, can be optimized globally by lifting the problem into a higher dimensional space where it admits a convex representation. Recent works, e.g. [41, 40], discuss
how to discretize the continuous formulation accurately.
We therefore expect to be able to solve all nonconvex majorizers Euk that are a discretization of (10) to (near)-global optimality, allowing us to consider highly non-trivial majorizing
functions. Similar relaxation methods exist in the discrete community via graph cuts for
Markov random fields, see [32, 11, 28] and the references therein.

1.2

Organization of this work

This work introduces an optimization algorithm for the sum of a function and a composite
function, which iteratively minimizes a nonconvex majorizing function (3). The algorithm
is detailed and discussed in Section 2 and basic properties are discussed in the first part
of Section 3. The second part of Section 3 then extends these basic properties to a global
convergence under the KL-property and uniqueness of R-minimizing solutions. Several generalizations and implementation details follow in Section 4. Finally, Section 5 shows some
promising numerical results on synthetic examples where the proposed algorithm is able to
find better minima than competing first-order methods, while being much more efficient than
methods from global optimization applied to the discussed problem class (2). We then close
Section 5 with an application to depth super resolution from noisy time-of-flight data.

2

The General Principle

Before we begin the formal introduction of the necessary context and provide convergence
and basic properties in their full generality it is instructive to reduce the problem formulation
to a very simple test case.
Let us consider the standard Jacobi-iteration:
uk+1 = D−1 (f − (A − D)uk )

(11)

which solves the linear equation Au = f for symmetric A ∈ Rn×n and its diagonal D by
successively minimizing the function
1
1
Euk (u) = hu, Du + (A − D)uk − f i − huk , (A − D)uk i,
2
2

(12)

which is a majorizer to E(u) = 12 hu, Au − f i, if D − A is positive definite. Now we would like
to solve the nonlinear equation system Aρ(u) = f for some function ρ : Rm → Rn . And we
do the same as before and apply our previous majorizer to ρ(u):
1
1
uk+1 = arg minhρ(u), Dρ(u) + (A − D)ρ(uk ) − f + hρ(uk ), (A − D)ρ(uk )ii.
2
2
u

(13)

If ρ is separable, then these problems can still be solved efficiently in each dimension, thereby
iteratively solving Aρ(u) = f . While this illustrates the main idea of our algorithm, the
situation becomes even more interesting if an additional regularization R makes a substitution
like z = ρ(u) impossible.
7

2.1

The Algorithm

Now we are ready to formulate the algorithm in full generality.
We consider the task of minimizing functions E : Rn → R ∪ {+∞} =: R and define the
domain of E by dom E =: {u ∈ Rn | E(u) < ∞}. We assume our energy to be proper,
that is E(u) > −∞ ∀u ∈ Rn , and its domain to be nonempty. We call a function lowersemicontinuous if lim inf u→ū E(u) ≥ E(ū).
We consider an optimization problem of the form
min E(u) = G(ρ(u)) + R(u),

u∈Rn

(14)

with G : Rm → R being smooth and proper, ρ : Rn → Rm continuous and R : Rn → R
proper and lower semi-continuous.
We define the Bregman distance of two vectors u, v ∈ Rn relative to a smooth and convex
’kernel’ function h by
Dh (u, v) = h(u) − h(v) − h∇h(v), u − vi.
We choose a step size τ > 0 to be discussed later, and then apply the following iterative
scheme:
Main Algorithm:
1
uk+1 ∈ arg min Dh (ρ(u), ρ(uk )) + h∇G(ρ(uk )), ρ(u) − ρ(uk )i + G(ρ(uk )) + R(u) (15)
u∈Rn τ

The use of a Bregman distance is an immediate generalization of induced squared norms, e.g
via h(u) = 12 ||u||22 , which allows us a natural inclusion of logarithmic or quadric G or vectorvalued step sizes. The connection to induced squared norms follows directly. Considering
any kernel function h = 21 ||u||2A , for a symmetric matrix A, we see that
1
1
1
Dh (u, v) = hAu, ui − hAv, vi − hAv, u − vi = ||u − v||2A .
2
2
2
Returning to the Jacobi example from before, we now see in particular that setting G(v) =
− f i, R(u) = 0 and h(v) = 12 ||u||2D exactly recovers the nonlinear Jacobi updates.
In practice this algorithm is feasible even if (15) is only locally solvable. However it is
especially interesting if (15) can actually be solved globally. In our applications we mainly
consider three interesting cases for this, although our theoretical analysis in later chapters is
not necessarily limited to those.
First, if ρ and R are Lipschitz with a bounded domain and separable, in the sense that
n can be written as ρ(u) = (ρ (u ), . . . , ρ (u )) and R : Rn → R can be written
ρ : Rn → RP
1 1
n n
as R(u) = ni=1 ri (ui ), then
(15)
decomposes
into
one-dimensional
subproblems for each ui .
P
We use separable h(u) = m
h
(u
),
so
that
d
(u
,
v
)
=
h
(u
)
−
hi (vi ) − h0i (vi )(ui − vi )
i
i
i
i
i
i
h
i
i=1
and find that the majorizer decouples so that
1
2 hv, Av

∂G(ρ(uk ))
1
(ρi (ui ) − ρi (uki )) + ri (ui ).
uk+1
∈ arg min dhi (ρi (ui ), ρi (uki )) +
i
τ
∂u
ui
i
8

(16)

Figure 4: Proximity relative to ρ can be crucial during minimization. Initialization marked in
black, global minimum in dark blue. Our majorizer, i.e. (15) and its minimizer are marked in
green. A majorizer that measures proximity relative to u, i.e. (6), shown in yellow. The left
figure shows a single step, the right figure shows the algorithm output and final majorizers.
These univariate nonconvex problems can be solved very efficiently and in parallel by uniform
grid searches or more elaborate exhaustive branch-and-bound strategies, due to the Lipschitz
properties R and ρ and the bounded domain.
A particularly interesting and practically relevant case are energies of the form


m
n
n
X
X
X
E(u) =
Fi 
ρij (uj ) +
ri (ui ),
(17)
i=1

j=1

i=1

where we have ρij : R → R and ri : R → R, Fi : Rm → R and we again assume a bounded
domain. These models appear naturally in several nonlinear regression tasks. But again, the
problem can be decomposed into one-dimensional
and we apply our algorithm,
Pn
P subproblems
v
F
(
as the subproblems decouple if we set G(v) = m
j=1 ij ) and ρ = (ρ11 , . . . , ρmn ).
i=1 i
Remarkably, both of the above examples still yield (near)-globally solvable subproblems, if
the separable regularization is replaced by a suitable penalty on the gradient of the unknown.
While such subproblems are nonconvex and non-separable they can still be solved efficiently
with the lifting techniques discussed in the context of equation (10). We detail these types
of problems in Section 4.3.

2.2

Special Cases

It is important to note several cases, where the method reduces to simpler approaches: First,
if ρ is the identity, then we immediately recover a non-composite problem, the setting of
[9]. If ρ is invertible, then we can minimize over z with the regularizer R(ρ−1 (z)) and again
recover a non-composite problem. Further, if G is separable as well, then it would be easier to
take the whole problem directly as a nonconvex majorizer, which would converge in a single
iteration. If the regularizer R is zero, then the algorithm works fine, yet we would like to
highlight that is possibly easier to solve the minimization over G(v) first under the constraint
of v ∈ Im ρ (for separable ρ this would be an especially easy constraint), and then optimize
Dh (ρ(u), v ∗ ).

9

2.3

Proximity relative to inner function

The subproblems in (15) critically contain the term Dh (ρ(u), ρ(uk )), which measures the
proximity of ρ(u) to ρ(uk ). This is in contrast to standard schemes, which classically consider
the distance of u to uk , Dh (u, uk ). This choice, although motivated naturally by the previous
example of a nonlinear Jacobi scheme may seem odd at first glance, especially in comparison
to previous methods for composite optimization, e.g. [19, 18, 10] or [49].
However this formulation has several advantages, whenever the subproblems can still be
solved efficiently. The main advantage is the leverage we gain. By updating relative to ρ,
we are able to directly apply smoothness properties and subsequent descent lemmas for G,
easily finding a global majorizer in each step. Furthermore, our step size can now be chosen
analytically independent of ρ and all new iterates are feasible in the sense that uk+1 ∈ dom ρ
and ρ(uk+1 ) ∈ dom h.
To make a more intuitive argument, we also note that penalizing the direct proximity
between u and uk of course limits the updates uk+1 to a neighborhood of uk . If we are able
to solve subproblems to global optimality, then limiting our updates in a local area seems
unnecessary. If we penalize the proximity in ρ, then we only stay in a local area relative to G,
which is necessary, as we linearized G. But otherwise we allow for arbitrarily large updates
as long as ρ(u) is similar to ρ(uk ), which is of no issue, as we solve our subproblems globally.
By this approach we hope to find interesting stationary points globally and not just locally
in a neighborhood around the starting vector.
Figure 4 visualizes this behavior in 1D. Given G, a smooth version of min(u2 , λ), ρ(u) =
sin(u) and R(u) = α|u| we majorize around the black mark using h(u) = u2 . We see that the
proximity relative to ρ is critical for reaching the global minimizer.
It is quite instructive to compute both update steps for a linear composition example, i.e.
E(u) = F (Au) for A ∈ Rm×n . One can check that the updates relative to ρij (uj ) =: aij uj
in (17) then correspond to a gradient descent with diagonal preconditioning, whereas the
update in u directly would correspond to standard gradient descent.

3

Algorithm Discussion and Convergence

In the following section we will analyze convergence properties of the proposed algorithm. We
will specify the smoothness assumptions we make on G, discuss well-posedness of subproblems
and give a descent lemma that guarantees subsequence convergence. We will then make
further assumptions on tameness of the function to prove global convergence.
To evaluate first-order descent with analytical step sizes we need to assume some notion of
bound on change in the gradient of a function. A recent generalization of Lipschitz continuity
of the gradient of G [5, 9] defines the following property:
Definition 3.1 (L-smooth adaptable). A proper, lower-semicontinuous function G : Rn → R
is called L-smooth adaptable if there exist L > 0 and a convex differentiable kernel function
h : Rn → R so that Lh − G is convex.
We measure proximity relative to this kernel function h through the Bregman distance it
induces,
Dh (u, v) =: h(u) − h(v) − h∇h(v), u − vi.

10

As a consequence of the L-smooth adaptability property we have the following descent
lemma for smooth G.
Lemma 1 (Descent Lemma, [5, Lemma 1]). If the function G is L-smooth adaptable to a
convex differentiable kernel function h, then
G(u) − G(v) − h∇G(v), u − vi ≤ LDh (u, v)

∀u, v ∈ Rn .

Proof. Df (u, v) ≥ 0 ∀u, v ∈ Rn if and only if f is convex. As such, DLh−G (u, v) ≥ 0 which
yields DG (u, v) ≤ LDh (u, v), due to the additivity of the Bregman distance.

3.1

Basic Properties

We define the subproblem energy in the following by
Euk (u) =

1
Dh (ρ(u), ρ(uk )) + G(ρ(uk )) + h∇G(ρ(uk )), ρ(u) − ρ(uk )i + R(u).
τ

(18)

For the first part of this section, we now make the following assumptions:
Assumptions A:
• E is proper, lower-semicontinuous and bounded from below
• E is coercive
• G is L-smooth adaptable for some convex kernel function h
• Every iteration is solved sufficiently, so that Euk (uk+1 ) ≤ Euk (uk ) = E(uk )
Under these assumptions we will discuss the validity of (18) as a majorizer for E and the wellposedness of the minimization of Euk (u). Note that in practice, we often gain coerciveness by
considering functions E with a bounded domain. The most important assumption here is the
smoothness assumption on G given by its L-smooth adaptability. The fourth assumption is
very general and holds, for example, already when each sub-problem is solved only by finite
sampling.
Lemma 2 (Majorization Property). Under the assumptions A and given some τ <
is a valid majorizer, i.e. it fulfills the properties

1
L,

Euk (u)

• Euk (u) ≥ E(u) ∀u ∈ Rn
• Euk (u) = E(uk ).
Proof. A quick computation shows that Euk (uk ) is equal to E(uk ), as the Bregman distance
Dh (x, y) is zero if x = y. Now applying Lemma 1 to ρ(u) and ρ(uk ) gives
G(ρ(u)) ≤ G(ρ(uk )) + h∇G(ρ(uk )), ρ(u) − ρ(uk )i + LDh (ρ(u), ρ(uk )).
On the other hand, we have, due to

1
τ

> L,

Euk (u) ≥ LDh (ρ(u), ρ(uk )) + G(ρ(uk )) + h∇G(ρ(uk )), ρ(u) − ρ(uk )i + R(u).
Combining both inequalities gives the desired result.
11

Now let us consider the set of minimizers of Euk ,(18):


k
n
Mτ (u ) = ū ∈ R | ū ∈ arg min Euk (u)

(19)

u

Lemma 3. Under the assumptions A, the set Mτ (uk ) is non-empty and compact for any
uk ∈ dom E if τ < L1 .
Proof. We already know that E is coercive. However, as Euk is a valid majorizer for τ < L1 ,
Euk (u) ≥ E(u), it is itself coercive. Furthermore E and hence Euk is bounded from below. As
a result Euk (u) is lower semi-continuous and proper with bounded level sets. [52, Theorem 1.9]
now guarantees that inf Euk is finite and that the set Mτ (uk ) is non-empty and compact.

3.2

Descent Lemma

As usual for majorization-minimization algorithms, we gain a monotone decrease in the objective function:
Lemma 4 (Descent Lemma). If the assumptions A hold, then the energy E is monotonically
decreasing for all iterates uk if τ < L1 is chosen as step size. The descent rate is
E(uk+1 ) − E(uk ) ≤ −

1 − τL
Dh (ρ(uk+1 ), ρ(uk ))).
τ

(20)

Proof. Using Lemma 1 for G we insert ρ(uk+1 ) and ρ(uk ) so that
G(ρ(uk+1 )) − G(ρ(uk )) ≤ h∇G(ρ(uk )), ρ(uk+1 ) − ρ(uk )i + LDh (ρ(uk+1 ), ρ(uk )).
Furthermore, because every iteration is solved sufficiently, we know that
1
Dh (ρ(uk+1 ), ρ(uk ))+h∇G(ρ(uk )), ρ(uk+1 )−ρ(uk )i+R(uk+1 )+G(ρ(uk )) ≤ G(ρ(uk ))+R(uk ).
τ
Adding both inequalities yields
1
G(ρ(uk+1 )) − G(ρ(uk )) + R(uk+1 ) − R(uk ) ≤ (L − )Dh (ρ(uk+1 ), ρ(uk )),
τ
which is the desired result. Due to the convexity of h, the right-hand side is always nonnegative if τ < L1 .
From this we gain convergence in function values, subsequence convergence and some
notion of convergence speed. For clarity we define the enveloping sequence z k = ρ(uk ).
Corollary 1 (Descent Lemma Corollary). Under the assumptions of A, we see that
∗
• the sequence of function values (E(uk ))∞
k=1 converges to a limit E ,

• limk→∞ Dh (z k+1 , z k ) = 0,
kl ∞
• there exist converging subsequences (uki )∞
i=1 and (z )l=1 ,

• the sequence min1≤k≤N Dh (z k+1 , z k ) converges to 0 with order O( N1 ).
12

Proof. Summing the descent lemma for N iterates uk and telescoping gives
α

N
X

Dh (z k+1 , z k ) ≤ E(u1 ) − E(uN ) ≤ C

k=1
L
k ∞
for α = 1−τ
τ . (E(u ))k=1 is a monotone decreasing sequence, that is bounded as ∞ < inf E <
E(uk ) ≤ E(u1 ) and thus converging. We gain the existence of converging subsequences uki
due to these bounds and the lower-semicontinuity and coercivity of E. The continuity of ρ
allows us to extend this to the existence of converging subsequences z kl .
Concerning the convergence rate, define a minimal proximity over all iterates

µN =: min Dh (z k+1 , z k ),
1≤k≤N

so that

N
X

Dh (z k+1 , z k ) ≤

k=1

E(u1 ) − E(uN )
α

implies
µN ≤

E(u1 ) − E(uN )
C
≤ .
αN
N

For later use we define the set of all accumulation points of the sequence uk , generated
by our algorithm for a given starting vector u0 as
accum(u0 ) = {u ∈ Rn | lim ukl = u for a subsequence ukl }.
l→∞

(21)

This set is non-empty as the sequence is bounded and closed as a set of limit points.

3.3

Convergence Properties

In this section we want to prove further statements of convergence. Up to now, we only
gave assumptions on G and on the relative minimization of the subproblems. For arbitrary
nonconvex ρ and R we can thus not expect a global convergence of the sequence of iterates,
uk , mostly because the accuracy up to which each subproblem is solved has not been specified
yet. However, even if the subproblems are solved exactly, we need to assume some algebraic
properties of E.
We will first give an appropriate optimality condition for our subproblems and specify the
algebraic notion of ’tameness’ discussed previously. Under these assumptions we will show
a global convergence of the sequence z k = ρ(uk ) and convergence to critical points. This is
a natural convergence result as the distance of successive iterates is only measured relative
to Dh (ρ(u), ρ(uk )). It can be a conscious modeling choice to allow several equivalent critical
points u∗ to be found by a single run of the algorithm. However we will also see that the
choice of regularizer R directly controls a global convergence in uk , if the subproblems are
solved to global optimality.
From now on, we consider a standard limiting subgradient:

13

Definition 3.2 ([52, 8.3]). A function E : Rn → R has the subgradient v ∈ Rn at a point
ū ∈ dom E, if
E(u) − E(ū) + hv, u − ūi
lim inf
≥0
u→ū,
|u − ū|
u6=ū

ˆ
and we write v ∈ ∂E(ū).
v is further an element of the limiting subgradient at ū if sequences
i
i ). We write
ˆ
exist so that u → ū, while also E(ui ) → E(ū) and v i → v for v i ∈ ∂E(u
v ∈ ∂E(ū) in this case.
Note that in our case the set ∂Euk (u) is non-empty due to Lemma 3. Under this definition
we can apply Rockafellar’s generalization of Fermat’s theorem.
Theorem 1 (Fermat’s theorem [52, 10.1]). If a proper function E : Rn → R has a local
minimum at ū ∈ dom E, then 0 ∈ ∂E(ū). If E = f + g, for a smooth function g and a proper
function f , then −∇g(ū) ∈ ∂f (ū).
From now on we assume ρ to be continuously differentiable and denote its Jacobian by
Jρ . This allows us to reformulate our problem as the search for critical points ū ∈ dom E
that fulfill
− Jρ (ū)∗ ∇G(ρ(ū)) ∈ ∂R(ū),
(22)
as a consequence of Theorem 1 and the chain rule for smooth functions. We further call the
set of all points ū that fulfill this condition critE.
However considering just the subdifferential of arbitrary functions leaves too many pathological cases for successful analysis of convergence properties [17]. As such we follow recent
literature on nonconvex optimization and consider functions that further satisfy the KurdykaLojasiewicz property:
Definition 3.3 (Nonsmooth Kurdyka-Lojasiewicz property [9]). For a proper and lowersemi-continuous function E : Rn → R, we define its local Kurdyka-Lojasiewicz property
(KL) at a point ū ∈ dom E by the attribution that there exist η > 0, ϕ ∈ C 0 [0, η) ∩ C 1 (0, η)
with ϕ(0) = 0, ϕ concave, ϕ0 > 0 and a neighborhood U (ū), so that
ϕ0 (E(u) − E(ū)) dist (0, ∂E(u)) (u) ≥ 1,
for all u ∈ U (ū) with E(ū) < E(u) < E(ū) + η.
If E is for example semi-algebraic, then it satisfies the KL-property at any ū ∈ dom E.
The graph of the limiting subdifferential of a proper semi-algebraic function E : Rn → R has
dimension n and thus E has a finite number of critical points [17]. We note that any function
definable in an o-minimal structure satisfies the KL-property [27]. Further, the property can
be uniformized to yield
Lemma 5 ([7, Lemma 6]). Let Ω be a compact set and consider a proper, lower-semicontinuous
function E : Rn → R. If E is constant on Ω and satisfies the KL-property at every point in
Ω, then there exist ε > 0, η > 0, ϕ ∈ C 0 [0, η) ∩ C 1 (0, η) with ϕ(0) = 0, ϕ concave, ϕ0 > 0
such that for all ū in Ω the uniformized KL-property,
ϕ0 (E(u) − E(ū)) dist (0, ∂E(u)) ≥ 1,
holds for all u ∈ Rn with dist(u, Ω) < ε and E(ū) < E(u) < E(ū) + η.
Now we are ready to extend our set of assumptions.
14

Assumptions B:
• E is continuous on dom E and meets the KL-property
• h is m-strongly convex
• ρ is a continuously differentiable map
•

1
τh

− G is L0 -smooth on any bounded subset

• Every iteration determines a uk+1 , so that 0 ∈ ∂Euk (uk+1 ).
These extended assumptions now allow us to prove the following statements:
1. The sequence of z k = ρ(uk ) converges globally to a value z ∗ .
2. All limit points of subsequences (uki )∞
i=1 are stationary points of E.
3. There is a correspondence between the limit point z ∗ and the accumulation points of
the iterates uk , given by z ∗ = ρ(u∗ ) for all u∗ ∈ accum(u0 ).
4. If the subproblems are solved globally, i.e. Euk (uk+1 ) ≤ Euk (u), then all accumulation
points are R-minimizing solutions to ρ(u) = z ∗ . If the R-minimizing solution to the
equation ρ(u) = z ∗ is unique, then the sequence uk converges globally.
We will see that the basic ingredients of the proof are basically as described in the recent
papers [2, 9], however special care has to be taken as all estimates of slope and value descent
are only relative to the ”enveloping” sequence of z k = ρ(uk ). Yet by appropriate arguments,
we can show that this is sufficient to prove our results.
We will see in the proof that the most demanding properties in assumptions
 B are used to
prove a bound on the slope of iterates, i.e. the inequality dist 0, ∂E(uk+1 ) ≤ c||z k+1 − z k ||
for some c > 0. If this inequality might be guaranteed otherwise, then global convergence still
follows. The assumptions on optimality and bounds on h and G reduce to the assumptions
in [9], for the non-composite case ρ = Id.
Lemma 6 (Slope bound by enveloping iterates). If the assumptions A and B hold, then
C < ∞ exists, so that
dist(0, ∂E(uk+1 )) ≤ C||z k+1 − z k ||.
(23)
Proof. Consider the optimality condition of the update equation, as all subproblems are
solved exactly:


1
1
k+1
k
k+1
k+1 ∗
k
0 ∈ ∂R(u ) + (Jρ (u )) ∇G(ρ(u )) + ∇h(ρ(u )) − ∇h(ρ(u ))
τ
τ
and reformulate to

k+1 ∗
(Jρ (u )) ∇G(ρ(uk+1 )) − ∇G(ρ(uk ))



1
1
k+1
k
− ∇h(ρ(u )) + ∇h(ρ(u )) ∈ ∂R(uk+1 ) + (Jρ (uk+1 ))∗ ∇G(ρ(uk+1 ) .
τ
τ
15

Now we see that the left hand side is an element of ∂E(uk+1 ), so that we can estimate its
norm by




1
1
k+1
k+1 ∗
k+1
dist(0, ∂E(u )) ≤ ||Jρ (u ) ||op ∇
h − G (z
)−∇
h − G (z k ) .
τ
τ
n
By Lemma 4, we know that (z k )∞
k=1 is a bounded subset of R . As such we can apply our
1
0
assumption on the L smoothness of τ h − G. Furthermore ||Jρ (u)∗ || is Lipschitz with some
constant M > 0 as C 1 implies Lipschitz continuity on all compact subsets, which concludes
the proof.

Before we now come to the main theorem, we first collect a few properties of the set
accum(u0 ), that will allow the application of Lemma 5.
Lemma 7. E is constant and finite on accum(u0 ), i.e. E(u) = E(v) < ∞ ∀u, v ∈ accum(u0 )
and we have
lim dist(uk , accum(u0 )) = 0.
(24)
k→∞

kl
Proof. (uk )∞
k=1 is bounded due to Corollary 1. We choose a subsequence (u ) so that
liml→∞ ukl = u∗ . From the continuity of E on (ukl ) we infer liml→∞ E(ukl ) = E(u∗ ). We
further know from Corollary 1 that the sequence of function values itself is convergent to
some value E ∗ , so that E is finite and constant on all these limit points. (24) is true for all
bounded sequences.

The following proof is now only a slight adaptation of usual strategies for convergence
under the KL-property, for example in Attouch et al. [2] or [9, 45].
Theorem 2 P
(Global Convergence). Under the previous assumptions, the sequence z k has
k+1 − z k || < ∞, and converges to a limit z ∗ .
finite length, ∞
k=1 ||z
Proof. If the iterates meet z k+1 = z k ∀k > N for some index N , then the theorem would be
clear. If not, then we may choose an index l large enough so that both E(ul ) < E ∗ + η and
dist(ul , accum(u0 )) < ε are fulfilled. The positive constants ε and η are the ones required by
the KL-property of E w.r.t to the set accum(u0 ). On this set, E is constant and finite, as
discussed in Lemma 7. From Lemma 5, we then find that




ϕ0 E(uk ) − E ∗ dist 0, ∂E(uk ) ≥ 1
holds for any accumulation point u∗ ∈ accum(u0 ), as E(u∗ ) = E ∗ and for all uk with k > l.
Now we can apply Lemma 6 to find that


ϕ0 E(uk ) − E ∗ ≥

C||z k

1
.
− z k−1 ||

(25)

Further, we can consider the descent from Lemma 4 and apply that h is m-strongly convex
to obtain
E(uk ) − E(uk+1 ) ≥

1 − τL
m(1 − τ L) k+1
Dh (z k+1 , z k ) ≥
||z
− z k ||2 .
τ
2τ
16

(26)

Analogously to [2, 9], we use the concavity of ϕ to analyze the difference of function values
in ϕ:




∆k,k+1 =: ϕ E(uk ) − E ∗ − ϕ E(uk+1 ) − E ∗



≥ ϕ0 E(uk ) − E ∗ E(uk ) − E(uk+1 ) .
Inserting (25) and (26), we gain
∆k,k+1 ≥ C 0

||z k+1 − z k ||2
.
||z k − z k−1 ||

Now we are entirely in the setting of [9, Theorem 6.2] and likewise reformulate to
q
2 ||z k − z k−1 ||C 0 ∆k,k+1 ≥ 2||z k+1 − z k ||
√
and use the inequality 2 ab ≤ a + b to gain
2||z k+1 − z k || ≤ ||z k − z k−1 || + C 0 ∆k,k+1 .
Summing these inequalities for k = l + 1, . . . , n then yields
2

n
X

||z k+1 − z k || ≤

k=l+1

=
≤

n
X
k=l+1
n
X
k=l+1
n
X

||z k − z k−1 || + C 0

n
X

∆k,k+1

k=l+1

||z k+1 − z k || − ||z n+1 − z n || + ||z l+1 − z l || + C 0

n
X

∆k,k+1

k=l+1

||z k+1 − z k || + ||z l+1 − z l || + C 0 ∆l+1,n+1 ,

k=l+1

where the last inequality is gained by telescoping all ∆. Reinserting the definition of ∆l+1,n+1
and simplifying then results in
n
X

||z k+1 − z k || ≤ ||z l+1 − z l || + C 0 ϕ(z l − z l+1 ) − C 0 ϕ(z n − z n+1 ) < ∞

k=l+1

As ϕ is positive this implies that the whole sequence z k is a Cauchy sequence and converges
due to metric completeness.
Remark 1 (Strong convexity of h). The m-strong convexity might seem like a limiting
assumption, yet it is always possible to construct a kernel function h̃ that fulfills this property,
if G is L-smooth adaptable for some h (see also the related discussion in [9]). First if Lh − G
2
is convex, then L(h + w) − (G + Lw) is also convex for any function w, especially w = m
2 || · || .
Define h̃ = h + w, G̃ = G + Lw, and R̃ = R − w ◦ ρ. The new energy G̃ ◦ ρ + R̃ is equal to
the old formulation, but the pair (G̃, h̃) is convex with h̃ being m-strongly convex.
With this enveloping convergence result of Theorem 2 we can return to uk :
17

Corollary 2 (Convergence to critical points). All accumulation points u∗ ∈ accumE of
∗
(uk )∞
k=1 are stationary points of E, i.e. u ∈ critE and belong to the same enveloping sequence
z k so that z ∗ = ρ(u∗ ).
Proof. Combining the bound on the slope in Lemma 6 and global convergence of z k ’s from
Theorem 2 we immediately see that dist(0, ∂E(uk+1 )) → 0 as k → ∞. Furthermore, we know
that E(uk ) → E(u∗ ) = E ∗ . As such all subsequences uk fulfill the definition of the limiting
subdifferential and 0 ∈ ∂E(u∗ ). We know also that limk→∞ z k = limk→∞ ρ(uk ) = z ∗ . Let
u∗ ∈ accum(u0 ) be arbitrary with the sequence by ukl → u∗ . Due to continuity of ρ we have
ρ(u∗ ) = liml→∞ ρ(ukl ) = z ∗ .
This result shows the connection between the enveloping ”proper” sequence of gradient
steps z k , which converges globally, due to the KL-property and the sequence of actual update
steps uk . The algorithm converges to a stationary point of E and all accumulation points not
only have the same value in E, but also in G ◦ ρ, as G(ρ(u∗ )) = G(z ∗ ) ∀u∗ ∈ accum(u0 ). A
necessary consequence of this is that all accumulation points also have equal value in R, so in
summary they are elements of the set C = {u ∈ Rn | R(u∗ ) = R(u), ρ(u∗ ) = z ∗ }. Naturally,
if this set is a singleton, then the subsequence convergence extends to global convergence.
The cardinality of the set C is however difficult to check a-priori. Nevertheless it turns out
that if we finally also assume that the nonconvex subproblems are solved globally, then the
convergence result follows from the familiar notion of uniqueness of R-minimizing solutions,
defined for example in [55], in the following way:
Definition 3.4. A vector u∗ ∈ Rn is called R-minimizing with respect to a solution set
{u ∈ Rn | F (u) = v} of an operator F : Rn → Rm and a vector v ∈ Rm , if F (u∗ ) = v and
R(u∗ ) ∈ min {R(u) | u ∈ Rn , F (u) = v} .
Now the uniqueness of such a vector directly corresponds to global convergence if the
subproblems are solved globally:
Theorem 3 (Global Convergence). If the subproblems are solved to global optimality, i.e.
Euk (uk+1 ) ≤ Euk (u) ∀u ∈ Rn , then
u∗ ∈ accum(u0 ) ⇒ u∗ is R-minimizing w.r.t ρ(u) = z ∗ .

(27)

In particular, if the R-minimizing element w.r.t. ρ(u) = z ∗ is unique, then the sequence of
iterates uk converges globally.
Proof. Let u ∈ {u ∈ Rn | ρ(u) = z ∗ } and u∗ ∈ accum(u0 ) be arbitrary. Rewriting the
optimality assumption Euk (uk+1 ) ≤ Euk (u) ∀u ∈ Rn results in the inequality
1
1
R(uk+1 )−R(u)+ Dh (ρ(uk+1 ), ρ(uk ))− Dh (ρ(u), ρ(uk ))+h∇G(ρ(uk )), ρ(uk+1 )−ρ(u)i ≤ 0.
τ
τ
Taking the limit of l → ∞ for a subsequence ukl → u∗ with the knowledge that liml→∞ ρ(ukl ) =
z ∗ then reveals that R(u∗ ) ≤ R(u), showing that u∗ ∈ accum(u0 ) is an R-minimizing solution
to ρ(u) = z ∗ . If in particular, the set of R-minimizing solutions is already a singleton, then
the result follows, as the set accum(u0 ) nonempty due to Corollary 1.
18

Example 1. As an example, consider a simple periodic function ρ : Rn → Rn , ρ(u) =
(sin(u1 ), . . . , sin(un )) and R = || · ||p for p > 0. The R-minimizing solution to ρ(u) = z is
then unique for any z ∈ Im(ρ). To see this consider that sin(x) is bijective on [− π2 , π2 ]. For
any level set zi ∈ [−1, 1] we can find a unique element ui in this interval [− π2 , π2 ], so that
sin(ui ) = zi . Due to the strict monotonicity of || · ||p on either R+ or R− , any other element
ui that fulfills ρ(ui ) = zi must have a greater function value in R.

4

Implementation Details

This section will focus on several interesting special cases of our general composite model (2)
and also discuss possible choices for the kernel function h.

4.1

Modeling

For several implementation examples it will be convenient to be a bit more specific with our
choices of G, ρ and R. One example is the natural extension to several additive terms,


m
n
n
X
X
X


E(u) =
Gi
ρij (uj ) +
rj (uj ),
(28)
i=1

j=1

j=1

which was already mentioned briefly in (17). This formulation is interesting due to its
straightforward interpretation as a way to optimize a function with n measurements of linear combinations of our m variables. Hence the task relates to nonlinear regression models
and imaging with nonlinear measurements. It is also a natural discretization of a general
nonlinear integral operator as defined for example
in [51,
Pn 4].
P
n
m×n , defined
G
(
However, setting G : Rm×n → R, G(v) = m
i
j=1 vij ) and ρ : R → R
i=1
component-wise by ρij (uj ), for univariate functions Gi , ρij , rj , we see that this is just an
instance of the general composite model2 . The maximal generalization would be achieved by
taking ρij : Rq → Rp , Gi : Rp → R, rj : Rq → R, although in practice q would have to be
quite small if we wanted to solve the subproblems by an exhaustive search.
PmWriting out the majorizer to (28) with univariate functions under the assumption that
i=1 Li hi − Gi is a convex function (as required for L-smooth adaptability, Definition 3.1)
gives


n
n
m
m
n
X
X
X
X
X
k
0
k 
Euk (u) =
Dhi (ρij (uj ), ρij (u )) +
Gi
ρij (uj )
ρij (uj ) +
rj (uj ), (29)
i=1

i=1

j=1

j=1

j=1

up to constant terms. This reveals
Pthat the majorization function is separable if each hi is
chosen separable so that hi (u) = nj=1 hij (uj ), as the Bregman distance to these hi is then
also separable and the summation over all m parameters can be exchanged with the summation over all n ’measurements’. The resulting m independent 1D dimensional subproblems
can then be solved efficiently.
2

Interestingly the univariate case of (28) exactly recovers functions described in the Kolmogorov–Arnold
representation theorem, although one would of course neither expect such a representation to be efficiently
computable nor the outer functions to be L-smooth adaptable.

19

Remark 2. This generalization is not only interesting for regression-type problems, where
the outer sum naturally sums over all samples and the inner sum over a superposition of
parametrized functions, but also for any sort of problem where it would make sense to split
a function into the composition of a function and a super-position of simpler functions. As
an example consider the 1-dimensional polynomial problem
P (x) =

p
X

!2
i

ai x − f

+

q
X

bi xi .

(30)

i=0

i=0

P
While it would be natural to choose ρ : R → R, ρ(x) = pi=0 ai xi −f , i.e the inner polynomial,
another possibility would
be to choose ρ : R → Rp , ρ(u) = (a0 , . . . , an xn ) and likewise to set
P
p
G : Rp → R, G(v) = ( j=1 vj − f )2 . A separable majorizer for this G using (29) would lead
to different subproblems than before.
An interesting fact about the general composite model (2) is that we are actually allowed
a great deal of freedom, as both G and ρ can be nonconvex. It is possible to insert any
invertible function f and its inverse f −1 on dom ρ and solve the equivalent problem with
G̃ = G ◦ f and ρ̃ = f −1 ◦ ρ. As an example, consider the following model


n
Y
E(u) = F 
gj (uj ) ,
j=1

where we have
of parametrized functions gj : R → R+ . We can set G : Rn →
 a product
P
R, G(v) = F exp( nj=1 vj ) and ρj : R → R, ρj (uj ) = log(gj (uj )), and recover the additive
superposition of parameters in (28).
We may freely use these possibilities due to Corollary 2. The proposed algorithm converges
to the set of stationary solutions of E. This result is independent of the actual decomposition
of E into G, ρ, R as long as the chosen triple G, ρ, R still fulfills all required conditions.

4.2

Computation of Bregman distance

Up to now, we always considered an abstract pair (G, h) fulfilling the conditions that both
functions are smooth and Lh − G and h are convex. Now we will detail several tangible
instances of these functions.
The trivial case is present when G is a concave function. We are then allowed to choose
an arbitrary convex function h, as −G is itself convex. A natural choice is to choose h as a
linear function, as its Bregman distance then vanishes,
Dh (u, v) = hh, ui − hh, vi − hh, u − vi = 0.
The resulting majorizer,
Euk (u) = h∇G(ρ(uk )), ρ(u) − ρ(uk )i + G(ρ(uk )) + R(u),

(31)

is an instance of iterative reweighting as discussed in [48, 47]. If R is convex and ρ is
coordinate-wise convex, then the majorizer is even convex. When R is a convex function and
20

ρ is an affine function, then we recover an instance of the difference of convex functions (DC)
algorithm [58].
The standard case is present when G is L-smooth. We then choose h = 21 || · ||2 and
recover the usual Euclidean distance measure via Dh (u, v) = 12 ||u − v||2 . Note that G can be
L-smooth without being convex, for example when considering a smooth truncated quadratic
function [1].
However, even when G is L-smooth, more advantageous kernel functions h might exist.
Consider the function G : Rn → R,
1
G(v) = ||Av − f ||2 ,
2

(32)

for a matrix A ∈ Rm×n . The function is L-smooth with L = ||AT A||op . However we can also
inspect Lh − G via its second derivative condition3 ,
L∇2 h(v) − AT A  0 ∀v ∈ Rn .

(33)

We could of course choose h = G, as (32) is convex, thereby solving the original problem
in each subproblem, but we are looking for functions h so that the subproblems are easy to
solve. Such a choice is presented by h = 12 || · ||2D with a diagonal matrix D. Choosing D so
that D − AT A  0 yields a diagonal preconditioning - we intrinsically find vector-valued step
sizes by an appropriate choice of h.
An important and motivating property of the L-smooth adaptable property is however
the inclusion of logarithmic functions, most prominently the Kullback-Leibler divergence as
possible terms for G, even though this function is not L-smooth at all [5]. Defining the
function G : Rn → R:


m
X
f
G(v) =
(Au)i − fi + fi log
,
(34)
(Au)i
i=1

Rm×n
+

forPA ∈
an appropriate kernel function on Rn++ is given by h(v) =
and f ∈
n
− i=1 log(vi ). [5, Lemma 7] reveals that the appropriate constant is L = ||f ||1 .
A general observation, see [49, Example 33], is that once we have gained a Bregman
distance Dh from h, we may actually use a whole family of kernel functions hk as long as
they majorize h while being convex,
Rn++

hk ∈ {hk ∈ C 1 (Rn ) | hk − h convex }.
The induced Bregman distance then fulfills
Dhk (u, v) ≥ Dh (u, v) ∀u, v ∈ Rn .
A specific instance of this observation is choosing h first and then implementing adaptive
step-sizes in this fashion by varying hk .
This is of course only a short overview of possible pairs (G, h), further examples can be
found in [5, 9, 49].
3

We follow the notation that A is positive semi-definite if A  0.

21

4.3

Models related to total variation regularization

Models that include the total variation norm as a regularization are ubiquitous in imaging
tasks [12, 53]. As such we will devote this section to these special instances, which have been
a major motivation for our work.
First, we now note that functions containing total variation terms,
E(u) =

n
X

gi (ui ) + ||Du||1 ,

(35)

i=1

where D denotes a finite-difference gradient and gi are continuous functions are a special case
of the following form
n
X
E(u) =
νi (ui , (Du)i ),
(36)
i=1

Rd

with continuous functions ν : R ×
→ R that are convex in their second argument. These
functions are natural discrete representations of the continuous model (10), which can be
efficiently solved by functional lifting [50]. But as soon as the function ν is not separable into
ν1 , . . . , νn , lifting approaches are not applicable anymore.
However we would especially like to apply lifting approaches to functions of the form
E(u) = G(ρ(u)) +

n
X

γi ((Du)i ),

(37)

i=1

with ρ separable as before, G L-smooth adaptable and γi convex. While the lifting scheme of
[50] is not applicable due to the non-separability of G, this is nevertheless a special instance
of our
Pgeneral problem (2). Indeed we can write down the majorizer, assuming a separable
h = ni=1 hi can be found, as
Euk (u) =

n
X

k

k

k

dhi (ρ(u), ρ(u )) + h∇G(ρ(u )), ρ(u )i +

i=1

n
X

γi ((Du)i ),

(38)

i=1

up to constants. Now this majorizer is in turn a particular instance of (36) and we can now
solve (37) by iteratively solving the lifting subproblems.
Furthermore, we can even exchange convexity of γi for differentiability. For arbitrary γi
that are L-smooth adaptable, we can linearize the second term as well, in full analogy to
(29), giving the majorizer
Euk (u) =

n
X

dhG (ρ(u), ρ(uk )) + h∇G(ρ(uk )), ρ(uk )i

i=1

+

n
X

dhγi ((Du)i , (Duk )i ) + γi0 ((Duk )i )(Du)i .

i=1

which is again an instance of (36). For concave γi this is particularly attractive as we can
choose hγi a as linear functions and just keep the linearization in full analogy to iterative
reweighting as discussed in the previous subsection in (31).
22

While this approach greatly increases the applicability of lifting schemes, it is important
to keep in mind that previous global optimality considerations for lifting schemes [50] do
not translate to these generalized problems. From Section 3 we only gain convergence to
stationary points. We will use the next section to analyze the quality of solutions that we
receive numerically.

4.4

Inertia

A small side note to the previous investigations that is nevertheless quite interesting in the
context of nonconvex optimization is inertia. Once we have (15), we can just as well consider
1
Euk ,uk−1 (u) = Dh (u, uk ) + h∇G(ρ(uk )), ρ(u) − ρ(uk )i + G(ρ(uk ))
τ

β
+ R(u) +
Dh (ρ(u), ρ(uk )) − Dh (ρ(u), ρ(uk−1 )) ,
τ

(39)

inserting an inertial term into the generalized forward-backward equation analogous to [46].
Inertia can be quite valuable for first-order optimization, especially as we allow G to be
nonconvex, but only utilize its gradient, i.e local information in each step. In practice we
observe that inertia can sometimes help the algorithm to reach better minima or speed up
the initial convergence for badly conditioned G. Also spurious stationary points can often
be overcome. Furthermore, the additional cost of solving (39) is minuscule compared to the
non-inertial variant.
However, due to the non-existence of the triangular inequality for Bregman distances, we
cannot bound the iterations by a Lyapunov function in general as in previous work [46] and
continue the proof of convergence with this Lyapunov function as a majorizer analogous to
Section 3. For the special cases of induced squared norms, i.e. h = ||u||2A , convergence still
follows by adapting Section 3.3 to the results of the recent work [45], but we omit a further
discussion.
In practice this modification still works well in many cases, setting β < 0.5. It is also
possible to backtrack in case of violations of Lyapunov function bounds.

5

Experimental results

In this section we analyze the proposed algorithm numerically. We will first consider a
synthetic example, where we will be able to compare the algorithm to other methods easily.
We will then move to an imaging application, the depth super-resolution from raw time-offlight data.

5.1

Synthetic experiments

We analyze the following energy
min

u∈Rn ,ui ∈[a,b]

Ff (Aρ(u)) + R(u),

(40)

where we have a boundedPinterval [a, b], an L-smooth adaptable function Ff ◦ A : Rn → R
n
∗
and regularizer R(u) =
i=1 r(ui − ui ) with r : R → R. ρ is chosen separable so that
23

(a)

(b)

(c)

Figure 5: Example of nonlinearities used in synthetic experiment, (a) x2 − 10 cos(2πx) (Rasx2
trigin’s function [42]), (b) random spline function with 12 queries and (c) 1+x
2.
ρ(u) = (ρ1 (u1 ), . . . , ρn (un )) with ρi : R → R, whom we will in general choose equal, and
omit the subscript. The nonlinearity r is aligned so that arg minx r(x) = 0. We first draw
u∗ ∈ [a, b]n at random and then set f = Aρ(u∗ ). We choose Ff as some measure of distance
between f and Aρ(u) and can thus guarantee that u∗ will be a global minimizer of (40).
Now we vary the difficulty of this possibly nonconvex optimization problem in two ways.
First we choose nonlinearities ρ, r as either
(1) Simple ρ(x) = exp(x), r(x) = x2
(2) Doable ρ(x) = x2 − 10 cos(2πx) [42], cf. Figure 5a, r(x) =

x2
,
1+x2

cf. Figure 5c

(3) Difficult ρ is a (coercive) piecewise cubic polynomial drawn by interpolating 12 values
in [a, b], r(x) = −sinc(x)
(4) Very Difficultρ is a (coercive) piecewise cubic polynomial drawn by interpolating 12
values in [a, b], r(x) = x2 − 10 cos(2πx).
This allows us to move from a nicely behaved, almost convex test case (1) to a nonconvex
problem with a well-behaved minimizer (2), adding further oscillations in (3) and finally
arriving at two very nonconvex functions in (4).
Then we vary the function G = Ff ◦ A. Note that this function critically determines the
interconnection of variables. If A is a diagonal matrix, then the problem is fully separable and
can by solved by n separate 1D optimizations with a single step of the algorithm, but if A is
a full matrix, then all variables are interdependent. Further, when A is a rectangular matrix,
then the system of nonlinear equations is under-determined and the function landscape is
(intuitively) not as well-behaved. Also, we are allowed to choose nonconvex functions F as
long as G is still L-smooth adaptable.
(a) Convex, local: Ff (v) = 21 ||v − f ||2 ,A ∈ Rn×n is a random matrix whose entries are
normally distributed relative to its diagonal. An appropriate
function is h(v) =
Pn kernel
1
2
T
j=1 |A A|ij
2 ||v||D , where D is a diagonal matrix with entries di =
P
(b) Convex, non L-smooth, local: Ff is the KL-divergence Ff (v) = ni=1 vi − fi log(vi ),
A ∈ Rn×n
P , is chosen as in (a). Here a kernel function is given by Burg’s entropy
h(v) = nj=1 − log(vj ) [5].

24

2e-13%

2.3e-11%

3.3e-14%

9.1e-14%

2 1.4e-13%

6.9%

1.6e-10%

2.5e-08%

3

2e-12%

4.5%

15%

4.1%

4

1.2%

3.5%

5.3%

3%

b

c

d

a

Grade of Nonlinearity

Grade of Nonlinearity

1

1 4.6e-14%

3.3e-11%

9.8e-15%

1.7e-14%

2 7.8e-14%

0.62%

2.3e-11%

3.3e-09%

3

2e-12%

4.5%

14%

3.6%

4

0.99%

3.4%

4.9%

2.9%

b

c

d

a

Difficulty in G

Difficulty in G

(a) Proposed method(15)

(b) Proposed (with inertia) (39)

Figure 6: Examination of different synthetic experiments. We increase the difficulty in G
from left to right and ∗in ρ, r from top to bottom, as detailed in Section 5. In each cell we
show the value E(ū)−E
, the value E(ū) reached by the algorithm minus the global minimum
Ẽ
∗
E , normalized by Ẽ, a (sampled) median of the energy values of the function, indicating the
quality of the solution relative to the overall energy landscape.
(c) Convex, full: Ff (v) = 21 ||v −f ||2 ,A ∈ Rm×n is a full random matrix with singular values
1
in [ log(n)
, 1]. m = n3 . Choose h as in (a).
(d) Nonconvex,
quadratic [1], i.e. a smoothed version of
P full: Ff is a smooth-truncated
2 , λ), A ∈ Rm×n is a full random matrix with singular values
min((v
−
f
)
Ff (v) = m
i
i
i=1
1
in [ log(n)
, 1], m = n3 . Choose h as in (a).
We run our method (15) without and with inertia, β = 0.4 (39). The subproblems in each
iterations are fully separable, so we solve the 1D problems in parallel by exhaustive search
with a sufficient amount of trial points and a parabolic refinement around the approximate
minimizer to desired precision. The refinement is a standard technique, e.g. [25] for 1D local
optimization and further references can be found, for example in book of Luenberger [38, pp.
217, 224]. This technique has also been used previously in imaging, for example, to refine
exhaustive search procedures in the context of quadratic decoupling for stereo in [33].
To mitigate the risk of lucky initializations, we run the algorithm 25 times with random
starting vectors for each test case and show the result which reached a median energy value.
We set n = 150 and [a, b] = [−3, 3], respectively [a, b] = [, 3] for the Poisson case and
implement the proposed method in MATLAB.
The results for all test cases can be found in Figure 6. We see that either increasing
the difficulty in G or the difficulty of the nonlinearity makes the overarching optimization
problem more difficult. Remarkably, our algorithm was able to find near-global optima for
25

10 5
10 5

10 0

10

10 0

-5

10
10

-5

-10

10 -10
10 -15

10

-20

10

0

10

1

10

2

10

3

(a) 1d

10 -15
10 0

10 1

10 2

(b) 3a

Figure 7: Convergence of various first-order methods for composite optimization to global
optimality. This figure shows problem 1d, i.e. a problem where R, ρ are easy and G difficult,
and 3a, the opposite case. Only our method provides favorable results for the difficult case
3a. Shown are the proposed method (15), proposed with inertia: (39), Adam: [31], gradient
descent: (1), forward-backward splitting: (4), Linear-Prox: (5), Prox-Linear: (6)
many test cases, especially the performance in row (2) is very good. However we see that the
increased oscillations in (3) eventually degrade the quality of solutions. Interestingly we also
notice that differences can appear even for convex functions G in (a) and (b). The squared
l2 norm in (a) seems to be easier to optimize globally, although the disparity to (b) is also
connected to the analytical step sizes, that we choose. Choosing larger stepsizes for (b) and
backtracking would recover a similar behavior to (a).
We now compare with other first-order nonlinear optimization methods that could possibly solve the same problem. Namely as mentioned in the related work section, gradient
descent (1), forward-backward splitting (4), the inner linearization, ’prox-linear’, (5) and the
outer linearization (6).
To fairly evaluate all majorizers we generally solve the subproblems in forward-backward
splitting (4), and outer linearization (6) to global optimality, again with exhaustive search
and parabolic fitting. For prox-linear (5), the subproblems do not decouple and we apply a
standard interior point solver in each iteration. We otherwise apply the same methodology
as before. We compare the convergence of all algorithms to the global minimum for two
characteristic cases, ’1d’ and ’3a’. While the first case denotes simple ρ, R and difficult G,
the second case denotes the opposite. We expect most methods to do well in the first case,
but the second case is not as clear. Figure 7a and Figure 7b show the results. It turns out
that indeed all methods can reliably solve the first case, Figure 7a, but only our method can
find near-optimal solutions for the second test case, Figure 7b. To compare our method to
a modern ’aggressive’ inertial variant that does not admit to a majorization-minimization
framework, we also include the Adam optimizer [31], however while this optimizer can find
better minima, it is still far off from the global solution in test case ’3a’.
We can go a step further and compare the proposed algorithm to global optimization
methods. See Figure 8 for a plot showing test case ’2c’ and the energetic difference to
the global minimum versus the runtime of each the algorithms, with time in a log-scale.
Previously mentioned algorithms are shown, as well as the MATLAB default implementations

26

10 5

Pattern Search
Genetic Algorithm
10 4

Prox-Linear

Algorithm Time

Multi Start
Simulated Annealing

10 3

Linear-Prox
10 2

Particle swarm
Proposed method

10 1

Gradient Descent
Forward-Backward Splitting
Adam

10 0
0

10

20

30

40

50

60

70

Distance to global optimum

Figure 8: Evaluation of various optimization methods for test case ’2c’, in terms of energetic
difference to global optimum vs. runtime.
of a genetic algorithm [22], particle swarm [30], pattern search [3], simulated annealing [26]
and multi-start methods [59]. All of these methods can reliably find global near-optimal
solutions for low dimensions, however in our setting of n = 150 these methods fail to find
a global minimizer within reasonable time constraints, as their efficiency decreases with the
number of variables. In contrast, our method exploits the structure of the objective function,
linearizing the convex outer function and solving the separable subproblems globally, and
scales well into higher dimensions.

5.2

Time-of-Flight Depth Reconstruction

Time-of-Flight cameras are used to recover depth images of a scene. They illuminate the
scene with a continuous wave and measure the time of flight of reflecting waves. A modern
hardware for this task are correlation photo-sensors, which directly measure the correlation
of the incoming wave with a reference wave [34]. The inversion of this correlation, however, is
highly non-linear and in practice often solved by assuming the incoming waves to be purely
sinosoidal, thereby including a well-known systematic error into the depth measurements,
especially at lower frequencies [37]. The inclusion of noise and a linear operator, e.g. here
a downsampling operator, into the reconstruction further inhibits the applicability of this
noise-free analytical inversion. A recent work on time-of-flight super-resolution [62] shows a
unified variational model. Further they show that modeling the incoming wave as arbitrary
periodic function and calibrating it thoroughly can increase the fidelity of the system. In
[62], the resulting nonconvex energy model is solved by alternating local optimization in
all variables. We will show that the problem can be solved with our approach and test on

27

-0.7

-0.8

-0.9

-1

-1.1

-1.2
-5

(a)

0

5

10

(b)

Figure 9: Ground Truth depth data [54] shown to the left and synthetic autocorrelation
function generated from trapezoidal signal to the right.

Figure 10: The four difference measurements yij , generated by equation (42) with subsampling and Gaussian noise.
synthetic data.
In the following, we will assume the following imaging model of a time-of-flight system


4πfi
2πj
ỹij = ai gi
u+
+ bi =: kij (u) + bi ,
(41)
λ
n
for measurement j in frequency fi and gi the 2π-periodic autocorrelation in frequency i that
is either calibrated or otherwise known, e.g as a cosine. ai is the amplitude in frequency fi ,
n the number of measurements in each frequency fi and λ the speed of light. To remove
the background illumination bi it is customary in Time-of-Flight literature to consider the
difference measurements
yij = kij (u) − ki,j+ n2 (u) =: ρij (u).
(42)
The problem of recovering a high-resolution depth image u from measurements yij can now
be stated as the energy minimization of
E(u) =

m−1
X n/2−1
X
i=0

||yij − Kρij (u)||2 + α||∇u||,

(43)

j=0

see also [62]. K is the imaging operator, which is here a downsampling operator. The total
variation regularization encourages a piecewise-constant depth solution.
The energy can be solved with the proposed method, because we can identify (43) with
the previously introduced special case of a sum of several composite terms (28) and a total
28

(a)

(b)

Figure 11: Classical closed form solution to depth recovery [34]. 90 MHz data to the left and
120 MHz to the right.
variation regularization (38). As such we can find a majorizer in each iteration that can be
solved by functional lifting.
In practice we solve all subproblems with sub-label accurate lifting as described in [41].
We initialize the algorithm with a constant depth of 1m and then iteratively update the
nonconvex majorizer and solve the sublabel-accurate lifting problem. We use a primal-dual
algorithm [14] to solve our subproblems and ’warm start’ each inner iteration with the primaldual variables from the previous step. We note that this relaxation approach can possibly
produce solutions that are convex combinations of global minimizers. As such we monitor
the energy of our inner iterations and terminate the algorithm if lifting cannot successfully
minimize the majorizer, i.e. if any iterates uk+1 were to violate Euk (uk+1 ) ≤ E(uk ), which
was postulated in Assumption A.
To verify this experimentally we generate synthetic data from a depth image of the Middlebury dataset [54] by applying (42). As a model for gi we use the autocorrelation of a
trapezoid signal. The autocorrelation signal and the ground truth depth are shown in Figure 9. We then apply downsampling by a factor of 2 to model the limited sensor size and
add significant Gaussian noise to model the sensitivity of common ToF sensors. We generate
two difference measurements in two frequencies, 90 MHz and 120 MHz. The ground truth
data covers a depth ranging from 0.5 to 6m. The resulting measurements are outside the
unambiguous range of both frequencies, so we expect a wrapping of data, which we want to
resolve using both frequencies. We visualize the resulting four difference measurements in
Figure 10. The noise level and severe data wrapping are apparent.
A classical inversion of the given data by the nonlinear closed form solution for sinusoidal
data [34] is shown for each frequency in Figure 11. The solution is however contaminated
by the nonlinear effects of noise and severe wrapping, note that we adjusted the colors for
visualization purposes. Further heuristics would be required in a next step to combine both
solutions to a final result, but we omit these due to the already significant distortion.
In contrast the solution by our algorithm is shown in Figure 12a. For reference, the solution to the algorithm, initialized with the ground truth data is also visualized in Figure 12b.
We see that the recovered solution is near-optimal. The algorithm can accurately unwrap and

29

(a)

(b)

Figure 12: Solution by the proposed algorithm with lifted subproblems, upsampling factor
of 2 to the left. To the right, proposed algorithm initialized with the ground truth, Figure 9,
for reference.
upsample most of the depth data and only small areas around the left chair are misidentified.
Remark 3. The presented model assumes the knowledge of the signal amplitude ai in each
frequency by some preceding algorithm to streamline the presentation. If this information
cannot be obtained robustly in practice, then the problem still falls the problem category
discussed in this paper, only the optimization variable v = [u, a] is then vector-valued in
depth and amplitude at each pixel (see the maximal generalization discussion in Section 4.1).
Yet, vector-valued variables can still be accounted for efficiently via the vectorial lifting shown
in the works [56, 35]. The overall algorithm remains unchanged, only the subproblems are
solved by vector-valued lifting instead.
We close this section by mentioning briefly that the presented composition of a matrix
and a nonlinear wrapping operator is not entirely unique to Time-of-Flight reconstruction.
A very related energy is present in nonlinear MRI reconstruction, see [60] for further reading.

6

Conclusions

In conclusion we proposed an optimization strategy for composite problems with simple,
but highly-nonlinear, inner functions and L-smooth adaptable outer functions. We construct
nonconvex majorizing functions and show that an iterative minimization of these functions
leads to the convergence of energy values under weak assumptions as well as the convergence of
the iterates to critical points of the energy under more restrictive assumptions. Our approach
has several attractive properties. It generates a set of feasible iterates and it is very easy to use
large step-sizes analytically, as these are independent of the Lipschitz properties of the inner
function. Our convergence results naturally extend previous work. In practice, we extensively
analyze the algorithm on synthetic examples, where the sub-problems can be solved globally
and show that it can find better minima than related methods. Lastly we show an intended
application. The use of recent functional lifting techniques to solve the nonconvex majorizer,

30

critically allows us to find visually appealing solutions to the complicated composite problem
of time-of-flight reconstruction from noisy low resolution data.

Acknowledgements
The authors acknowledge the support of the German Research Foundation (DFG) via the
research training group GRK 1564 Imaging New Modalities and the project ”Functional
Lifting 2.0”, as well as the support of the Daimler and Benz Foundation.

References
[1] M. Artina, M. Fornasier, and F. Solombrino, Linearly Constrained Nonsmooth
and Nonconvex Minimization, SIAM J. Optim., 23 (2013), pp. 1904–1937, https://
doi.org/10.1137/120869079.
[2] H. Attouch, J. Bolte, and B. F. Svaiter, Convergence of descent methods for
semi-algebraic and tame problems: Proximal algorithms, forward–backward splitting, and
regularized Gauss–Seidel methods, Math. Program., 137 (2013), pp. 91–129, https://
doi.org/10.1007/s10107-011-0484-9.
[3] C. Audet and J. Dennis, Analysis of Generalized Pattern Searches, SIAM J. Optim.,
13 (2002), pp. 889–903, https://doi.org/10.1137/S1052623400378742.
[4] C. Bardaro, J. Musielak, and G. Vinti, Nonlinear Integral Operators and Applications, Walter de Gruyter, Jan. 2003.
[5] H. H. Bauschke, J. Bolte, and M. Teboulle, A Descent Lemma Beyond Lipschitz
Gradient Continuity: First-Order Methods Revisited and Applications, Mathematics of
Operations Research, 42 (2017), pp. 330–348, https://doi.org/10.1287/moor.2016.
0817.
[6] A. Beck and M. Teboulle, A Fast Iterative Shrinkage-Thresholding Algorithm for
Linear Inverse Problems, SIAM Journal on Imaging Sciences, 2 (2009), pp. 183–202,
https://doi.org/10.1137/080716542.
[7] J. Bolte, S. Sabach, and M. Teboulle, Proximal alternating linearized minimization for nonconvex and nonsmooth problems, Math. Program., 146 (2014), pp. 459–494,
https://doi.org/10.1007/s10107-013-0701-9.
[8] J. Bolte, S. Sabach, and M. Teboulle, Nonconvex Lagrangian-Based Optimization:
Monitoring Schemes and Global Convergence, arXiv:1801.03013 [math], (2018), https:
//arxiv.org/abs/1801.03013.
[9] J. Bolte, S. Sabach, M. Teboulle, and Y. Vaisbourd, First Order Methods beyond
Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse
Problems, arXiv:1706.06461 [cs, math], (2017), https://arxiv.org/abs/1706.06461.
[10] S. Bonettini, I. Loris, F. Porta, and M. Prato, Variable Metric Inexact LineSearch-Based Methods for Nonsmooth Optimization, SIAM Journal on Optimization, 26
(2016), pp. 891–921, https://doi.org/10.1137/15M1019325.
31

[11] Y. Boykov, O. Veksler, and R. Zabih, Fast approximate energy minimization via
graph cuts, IEEE Transactions on Pattern Analysis and Machine Intelligence, 23 (2001),
pp. 1222–1239, https://doi.org/10.1109/34.969114.
[12] M. Burger and S. Osher, A Guide to the TV zoo, in PDE Based Reconstruction
Methods in Imaging, no. 2090 in Lecture Notes in Mathematics, Springer International
Publishing, Switzerland, 1 ed., 2013.
[13] A. Chambolle, D. Cremers, and T. Pock, A Convex Approach to Minimal Partitions, SIAM Journal on Imaging Sciences, 5 (2012), pp. 1113–1158, https://doi.org/
10.1137/110856733.
[14] A. Chambolle and T. Pock, A First-Order Primal-Dual Algorithm for Convex
Problems with Applications to Imaging, J Math Imaging Vis, 40 (2011), pp. 120–145,
https://doi.org/10.1007/s10851-010-0251-1.
[15] T. Chan, S. Esedoglu, and M. Nikolova, Algorithms for Finding Global Minimizers of Image Segmentation and Denoising Models, SIAM J. Appl. Math., 66 (2006),
pp. 1632–1648, https://doi.org/10.1137/040615286.
[16] G. Chen and R. Rockafellar, Convergence Rates in Forward–Backward
Splitting, SIAM J. Optim., 7 (1997), pp. 421–444, https://doi.org/10.1137/
S1052623495290179.
[17] D. Drusvyatskiy, Slope and Geometry in Variational Mathematics, PhD thesis, Cornell University, 2013.
[18] D. Drusvyatskiy, A. D. Ioffe, and A. S. Lewis, Nonsmooth optimization
using Taylor-like models: Error bounds, convergence, and termination criteria,
arXiv:1610.03446 [math], (2016), https://arxiv.org/abs/1610.03446.
[19] D. Drusvyatskiy and C. Paquette, Efficiency of minimizing compositions of convex
functions and smooth maps, arXiv:1605.00125 [math], (2016), https://arxiv.org/abs/
1605.00125.
[20] T. Frerix, T. Möllenhoff, M. Moeller, and D. Cremers, Proximal Backpropagation, arXiv:1706.04638 [cs], (2017), https://arxiv.org/abs/1706.04638.
[21] D. Gabay and B. Mercier, A dual algorithm for the solution of nonlinear variational
problems via finite element approximation, Computers & Mathematics with Applications, 2 (1976), pp. 17–40, https://doi.org/10.1016/0898-1221(76)90003-1.
[22] D. E. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning,
Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 1st ed., 1989.
[23] E. R. Hansen and G. W. Walster, Global Optimization Using Interval Analysis,
no. 264 in Monographs and textbooks in pure and applied mathematics, Marcel Dekker,
New York, 2nd ed., revised and expanded ed., 2004.
[24] D. R. Hunter and K. Lange, A Tutorial on MM Algorithms, The American Statistician, 58 (2004), pp. 30–37, https://doi.org/10.1198/0003130042836.
32

[25] W. Huyer and A. Neumaier, Global Optimization by Multilevel Coordinate Search,
Journal of Global Optimization, 14 (1999), pp. 331–355, https://doi.org/10.1023/A:
1008382309369.
[26] L. Ingber, Adaptive simulated annealing (ASA): Lessons learned, arXiv:cs/0001018,
(2000), https://arxiv.org/abs/cs/0001018.
[27] A. Ioffe, An Invitation to Tame Optimization, SIAM J. Optim., 19 (2009), pp. 1894–
1917, https://doi.org/10.1137/080722059.
[28] H. Ishikawa and D. Geiger, Segmentation by Grouping Junctions, in Proceedings of
the IEEE Computer Society Conference on Computer Vision and Pattern Recognition,
CVPR ’98, Washington, DC, USA, 1998, IEEE Computer Society, pp. 125–.
[29] R. B. Kearfott, Rigorous Global Search: Continuous Problems, Springer Science &
Business Media, Mar. 2013.
[30] J. Kennedy and R. C. Eberhardt, Particle Swarm Optimization, Proceedings of
the 1995 IEEE International Conference on Neural Networks (Conference proceedings),
4 (1995), pp. 1942–1948.
[31] D. P. Kingma and J. Ba, Adam: A Method for Stochastic Optimization,
arXiv:1412.6980 [cs], (2014), https://arxiv.org/abs/1412.6980.
[32] V. Kolmogorov and R. Zabin, What energy functions can be minimized via graph
cuts?, IEEE Transactions on Pattern Analysis and Machine Intelligence, 26 (2004),
pp. 147–159, https://doi.org/10.1109/TPAMI.2004.1262177.
[33] G. Kuschk and D. Cremers, Fast and Accurate Large-Scale Stereo Reconstruction
Using Variational Methods, in 2013 IEEE International Conference on Computer Vision
Workshops, Dec. 2013, pp. 700–707, https://doi.org/10.1109/ICCVW.2013.96.
[34] R. Lange, 3D Time-of-Flight Distance Measurement with Custom Solid-State Image
Sensors in CMOS/CCD-Technology, PhD thesis, University of Siegen, Siegen, June 2000.
[35] E. Laude, T. Möllenhoff, M. Moeller, J. Lellmann, and D. Cremers,
Sublabel-Accurate Convex Relaxation of Vectorial Multilabel Energies, in Computer Vision – ECCV 2016, Lecture Notes in Computer Science, Springer, Cham, Oct. 2016,
pp. 614–627, https://doi.org/10.1007/978-3-319-46448-0_37.
[36] A. S. Lewis and S. J. Wright, A proximal method for composite minimization, Math.
Program., 158 (2016), pp. 501–546, https://doi.org/10.1007/s10107-015-0943-9.
[37] M. Lindner and A. Kolb, Lateral and depth calibration of pmd-distance sensors, in
International Symposium on Visual Computing, Springer, 2006, pp. 524–533.
[38] D. G. Luenberger and Y. Ye, Linear and Nonlinear Programming, Springer, New
York, NY, 4th ed. 2016 edition ed., June 2015.

33

[39] J. Mairal, Optimization with First-order Surrogate Functions, in Proceedings of the
30th International Conference on International Conference on Machine Learning - Volume 28, ICML’13, Atlanta, GA, USA, 2013, JMLR.org, pp. III–783–III–791.
[40] T. Möllenhoff and D. Cremers, Sublabel-Accurate Discretization of Nonconvex
Free-Discontinuity Problems, Proceedings of the IEEE International Conference on Computer Vision, (2017), pp. 1183–1191.
[41] T. Möllenhoff, E. Laude, M. Moeller, J. Lellmann, and D. Cremers,
Sublabel-Accurate Relaxation of Nonconvex Energies, in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2016, pp. 3948–3956.
[42] H. Mühlenbein, M. Schomisch, and J. Born, The parallel genetic algorithm as
function optimizer, Parallel Computing, 17 (1991), pp. 619–632, https://doi.org/10.
1016/S0167-8191(05)80052-3.
[43] Y. Nesterov, Introductory Lectures on Convex Optimization, vol. 87 of Applied Optimization, Springer US, Boston, MA, 2004.
[44] Y. Nesterov, Gradient methods for minimizing composite functions, Math. Program.,
140 (2013), pp. 125–161, https://doi.org/10.1007/s10107-012-0629-5.
[45] P. Ochs, Unifying abstract inexact convergence theorems for descent methods and block
coordinate variable metric iPiano, arXiv:1602.07283 [math], (2016), https://arxiv.
org/abs/1602.07283.
[46] P. Ochs, Y. Chen, T. Brox, and T. Pock, iPiano: Inertial Proximal Algorithm
for Nonconvex Optimization, SIAM J. Imaging Sci., 7 (2014), pp. 1388–1419, https:
//doi.org/10.1137/130942954.
[47] P. Ochs, A. Dosovitskiy, T. Brox, and T. Pock, An Iterated L1 Algorithm for
Non-smooth Non-convex Optimization in Computer Vision, IEEE, June 2013, pp. 1759–
1766, https://doi.org/10.1109/CVPR.2013.230.
[48] P. Ochs, A. Dosovitskiy, T. Brox, and T. Pock, On iteratively reweighted algorithms for nonsmooth nonconvex optimization in computer vision, SIAM Journal
on Imaging Sciences, 8 (2015), pp. 331–372, https://doi.org/https://doi.org/10.
1137/140971518.
[49] P. Ochs, J. Fadili, and T. Brox, Non-smooth Non-convex Bregman Minimization:
Unification and new Algorithms, arXiv:1707.02278 [cs, math], (2017), https://arxiv.
org/abs/1707.02278.
[50] T. Pock, D. Cremers, H. Bischof, and A. Chambolle, Global Solutions of Variational Models with Convex Regularization, SIAM J. Imaging Sci., 3 (2010), pp. 1122–
1145, https://doi.org/10.1137/090757617.
[51] R. Precup, Methods in Nonlinear Integral Equations, Springer Science & Business Media, Mar. 2013.

34

[52] R. T. Rockafellar and R. J.-B. Wets, Variational Analysis, Springer Science &
Business Media, June 2009.
[53] L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal
algorithms, Physica D: Nonlinear Phenomena, 60 (1992), pp. 259–268, https://doi.
org/10.1016/0167-2789(92)90242-F.
[54] D. Scharstein, H. Hirschmüller, Y. Kitajima, G. Krathwohl, N. Nešić,
X. Wang, and P. Westling, High-Resolution Stereo Datasets with Subpixel-Accurate
Ground Truth, in Pattern Recognition, Lecture Notes in Computer Science, Springer,
Cham, Sept. 2014, pp. 31–42, https://doi.org/10.1007/978-3-319-11752-2_3.
[55] O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier, and F. Lenzen,
Variational Methods in Imaging, no. 167 in Applied Mathematical Sciences, Springer,
New York, 1 ed., 2009.
[56] E. Strekalovskiy, A. Chambolle, and D. Cremers, Convex Relaxation of Vectorial Problems with Coupled Regularization, SIAM Journal on Imaging Sciences, 7 (2014),
pp. 294–336, https://doi.org/10.1137/130908348.
[57] Y. Sun, P. Babu, and D. P. Palomar, Majorization-Minimization Algorithms in Signal Processing, Communications, and Machine Learning, IEEE Transactions on Signal
Processing, 65 (2017), pp. 794–816, https://doi.org/10.1109/TSP.2016.2601299.
[58] P. D. Tao and L. T. H. An, Convex analysis approach to dc programming: Theory,
algorithms and applications, Acta Mathematica Vietnamica, 22 (1997), pp. 289–355.
[59] Z. Ugray, L. Lasdon, J. Plummer, F. Glover, J. Kelly, and R. Martı́, Scatter
Search and Local NLP Solvers: A Multistart Framework for Global Optimization, INFORMS Journal on Computing, 19 (2007), pp. 328–340, https://doi.org/10.1287/
ijoc.1060.0175.
[60] T. Valkonen, A primal–dual hybrid gradient method for nonlinear operators with applications to MRI, Inverse Problems, 30 (2014), p. 055012, https://doi.org/10.1088/
0266-5611/30/5/055012.
[61] C. F. J. Wu, On the Convergence Properties of the EM Algorithm, The Annals of
Statistics, 11 (1983), pp. 95–103, https://doi.org/10.2307/2240463.
[62] L. Xiao, F. Heide, M. O’Toole, A. Kolb, M. B. Hullin, K. Kutulakos, and
W. Heidrich, Defocus deblurring and superresolution for time-of-flight depth cameras,
in 2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), IEEE,
2015, pp. 2376–2384.

35

