Peter Ochs

Jalal Fadili

Thomas Brox

Non-smooth Non-convex Bregman Minimization: Unification and new Algorithms

We propose a unifying algorithm for non-smooth non-convex optimization. The algorithm approximates the objective function by a convex model function and finds an approximate (Bregman) proximal point of the convex model. This approximate minimizer of the model function yields a descent direction, along which the next iterate is found. Complemented with an Armijo-like line search strategy, we obtain a flexible algorithm for which we prove (subsequential) convergence to a stationary point under weak assumptions on the growth of the model function error. Special instances of the algorithm with a Euclidean distance function are, for example, Gradient Descent, Forward-Backward Splitting, ProxDescent, without the common requirement of a "Lipschitz continuous gradient". In addition, we consider a broad class of Bregman distance functions (generated by Legendre functions), replacing the Euclidean distance. The algorithm has a wide range of applications including many linear and non-linear inverse problems in signal/image processing and machine learning.

Introduction

When minimizing a non-linear function on the Euclidean vector space, a fundamental strategy is to successively minimize approximations to the actual objective function. We refer to such an approximation as model (function). A common model example in smooth optimization is linearization (first order Taylor approximation) around the current iterate. However, in general, the minimization of a linear function does not provide a finite solution, unless, for instance, the domain is compact. Therefore, the model is usually complemented by a proximity measure, which favors a solution (the next iterate) close to the current iterate. For the Euclidean norm as proximity measure, computing the next iterate (minimizer of the sum of the model function and the Euclidean proximity measure) is equivalent to a Gradient Descent step, i.e. the next iterate is obtained by moving along the direction of the negative gradient at the current point for a certain step size.

Contributions and Related Work

Since sequential minimization of model functions does not require smoothness of the objective or the model function, non-smoothness is handled naturally. The crucial aspect is the "approximation quality" of the model function, which is controlled by a growth function, that describes the approximation error around the current iterate. Drusvyatskiy et al. [START_REF] Drusvyatskiy | Nonsmooth optimization using Taylor-like models: error bounds, convergence, and termination criteria[END_REF] refer to such model functions as Taylor-like models. The difference among algorithms lies in the properties of such a growth function, rather than the specific choice of a model function.

For the example of the Gradient Descent model function (linearization around the current iterate) for a continuously differentiable function, the value and the derivative of the growth function (approximation error) vanish at the current iterate. In this case, a line search strategy is required to determine a suitable step size that reduces the objective value. If the gradient of the objective function is additionally L-Lipschitz continuous, then the growth function satisfies a quadratic growth globally, and step sizes can be controlled analytically.

A large class of algorithms, which are widely popular in machine learning, statistics, computer vision, signal and image processing can be cast in the same framework. This includes algorithms such as Forward-Backward Splitting [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF] (Proximal Gradient Descent), ProxDescent [START_REF] Lewis | A proximal method for composite minimization[END_REF][START_REF] Drusvyatskiy | Error bounds, quadratic growth, and linear convergence of proximal methods[END_REF] (or proximal Gauss-Newton method), and many others. They all obey the same growth function as Gradient Descent. This allows for a unified analysis of all these algorithms, which is a key contribution of this paper. Moreover, we allow for a broad class of (iteration dependent) Bregman proximity functions (e.g., generated by common entropies such as Boltzmann-Shannon, Fermi-Dirac, and Burg's entropy), which leads to new algorithms. To be generic in the choice of the objective, the model, and the Bregman functions, the algorithm is complemented with an Armijo-like line search strategy. Subsequential convergence to a stationary point is established for different types of growth functions.

The above mentioned algorithms are ubiquitous in applications of machine learning, computer vision, image/signal processing, and statistics as is illustrated in Section 5 and in our numerical experiments in Section 6. Due to the unifying framework, the flexibility of these methods is considerably increased further.

Contributions and Related Work

For smooth functions, Taylor's approximation is unique. However, for non-smooth functions, there are only "Taylor-like" model functions [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF][START_REF] Noll | Convergence of non-smooth descent methods using the Kurdyka-Lojasiewicz inequality[END_REF][START_REF] Drusvyatskiy | Nonsmooth optimization using Taylor-like models: error bounds, convergence, and termination criteria[END_REF]. Each model function yields another algorithm. Some model functions [START_REF] Noll | A proximity control algorithm to minimize nonsmooth and nonconvex functions[END_REF][START_REF] Noll | Convergence of non-smooth descent methods using the Kurdyka-Lojasiewicz inequality[END_REF] could also be referred to as lower-Taylor-like models, as there is only a lower bound on the approximation quality of the model. Noll et al. [START_REF] Noll | Convergence of non-smooth descent methods using the Kurdyka-Lojasiewicz inequality[END_REF] addressed the problem by bundle methods based on cutting planes, which differs from our setup.

The goal of Drusvyatskiy et al. [START_REF] Drusvyatskiy | Nonsmooth optimization using Taylor-like models: error bounds, convergence, and termination criteria[END_REF] is to measure the proximity of an approximate solution of the model function to a stationary point of the original objective, i.e., a suitable stopping criterion for non-smooth objectives is sought. On the one hand, their model functions may be non-convex, unlike ours. On the other hand, their growth functions are more restrictive. Considering their abstract level, the convergence results may seem satisfactory. However, several assumptions that do not allow for a concrete implementation are required, such as a vanishing distance between successive iterates and convergence of the objective values along a generated convergent subsequence to the objective value of the limit point. This is in contrast to our framework.

We assume more structure of the subproblems: They are given as the sum of a model function and a Bregman proximity function. With this mild assumption on the structure and a suitable line-search procedure, the algorithm can be implemented and the convergence results apply. We present the first implementable algorithm in the abstract model function framework and prove subsequential convergence to a stationary point.

Our algorithm generalizes ProxDescent [START_REF] Drusvyatskiy | Error bounds, quadratic growth, and linear convergence of proximal methods[END_REF][START_REF] Lewis | A proximal method for composite minimization[END_REF] with convex subproblems, which is known for its broad applicability. We provide more flexibility by considering Bregman proximity functions, and our backtracking line-search need not solve the subproblems for each trial step.

The algorithm and convergence analysis is a far-reaching generalization of Bonettini et al. [START_REF] Bonettini | Variable metric inexact line-search based methods for nonsmooth optimization[END_REF], which is similar to the instantiation of our framework where the model function leads to Forward-Backward Splitting. The proximity measure of Bonettini et al. [START_REF] Bonettini | Variable metric inexact line-search based methods for nonsmooth optimization[END_REF] is assumed to satisfy a strong convexity assumption. Our proximity functions can be generated by a broad class of Legendre functions, which includes, for example, the non-strongly convex Burg's entropy [START_REF] Burg | The relationship between maximum entropy spectra and maximum likelihood spectra[END_REF][START_REF] Bauschke | Legendre functions and the method of random Bregman projections[END_REF] for the generation of the Bregman proximity function.

Preliminaries and Notations

Throughout the whole paper, we work in a Euclidean vector space R N of dimension N ∈ N equipped with the standard inner product •, • and associated norm | • |.

Variational analysis. We work with extended-valued functions

f : R N → R, R := R ∪ {±∞}. The domain of f is dom f := x ∈ R N | f (x) < +∞ and a function f is proper, if it is nowhere -∞ and dom f = ∅. It is lower semi-continuous (or closed), if lim inf x→x f (x) ≥ f (x) for any x ∈ R N . Let int Ω denote the interior of Ω ⊂ R N . We use the notation of f -attentive convergence x f → x ⇔ (x, f (x)) → (x, f (x)), and the notation k K → ∞ for some K ⊂ N to represent k → ∞ where k ∈ K.
As in [START_REF] Drusvyatskiy | Nonsmooth optimization using Taylor-like models: error bounds, convergence, and termination criteria[END_REF], we introduce the following concepts. For a closed function f : R N → R and a point x ∈ dom f , we define the slope of f at x by

|∇f |(x) := lim sup x→x, x =x [f (x) -f (x)] + |x -x| ,
where [s] + := max{s, 0}. It is the maximal instantaneous rate of decrease of f at x. For a differentiable function, it coincides with the norm of the gradient |∇f (x)|. Moreover, the limiting slope |∇f |(x) := lim inf

x f →x |∇f |(x) is key. For a convex function f , we have |∇f |(x) = inf v∈∂f (x) |v|, where ∂f (x) is the (convex) subdifferential ∂f (x) := v ∈ R N | ∀x : f (x) ≥ f (x) + x -x, v , whose domain is given by dom ∂f := x ∈ R N | ∂f (x) = ∅ . A point x is a stationary point of the function f , if |∇f |(x) = 0 holds. Obviously, if |∇f |(x) = 0, then |∇f |(x) = 0.
We define the set of (global) minimizers of a function f by Argmin

x∈R N f (x) := x ∈ R N | f (x) = inf x∈R N f (x) ,
and the (unique) minimizer of f by argmin

x∈R N f (x), if Argmin x∈R N f (x) consists of a single element.
As shorthand, we also use Argmin f and argmin f . Definition 1 (Growth function [START_REF] Drusvyatskiy | Nonsmooth optimization using Taylor-like models: error bounds, convergence, and termination criteria[END_REF]). A differentiable univariate function ω : R

+ → R + is called growth function if it satisfies ω(0) = ω + (0) = 0.
If, in addition, ω + (t) > 0 for t > 0 and equalities lim t 0 ω + (t) = lim t 0 ω(t)/ω + (t) = 0 hold, we say that ω is a proper growth function.

Concrete instances of growth functions will be generated for example by the concept of ψ-uniform continuity, which is a generalization of Lipschitz and Hölder continuity. Definition 2. A mapping F : R N → R M is called ψ-uniform continuous with respect to a continuous function ψ : R + → R + with ψ(0) = 0, if the following holds:

|F (x) -F (x)| ≤ ψ(|x -x|) for all x, x ∈ R N .
Example 3. Let F be ψ-uniform continuous. If, for some c > 0, we have ψ(s) = cs α with α ∈]0, 1], then F is Hölder continuous, which for α = 1 is the same as Lipschitz continuity.

In analogy to the case of Lipschitz continuity, we can state a generalized Descent Lemma: Lemma 4 (Generalized Descent Lemma). Let f : R N → R be continuously differentiable and let ∇f : R N → R N be ψ-uniform continuous. Then, the following holds

|f (x) -f (x) -∇f (x), x -x | ≤ 1 0 ϕ(s|x -x|) s d s for all x, x ∈ R N ,
where ϕ : R + → R + is given by ϕ(s) := sψ(s).

Proof. We follow the proof of the Descent Lemma for functions with Lipschitz gradient:

|f (x) -f (x) -∇f (x), x -x | = | 1 0 ∇f (x + s(x -x)) -∇f (x), x -x d s| ≤ 1 0 |∇f (x + s(x -x)) -∇f (x)||x -x| d s ≤ 1 0 ψ(s|x -x|)|x -x| d s = 1 0 ϕ(s|x -x|) s d s . Example 5. The function ω(t) = 1 0 ϕ(st)
s d s is an example for a growth function. Obviously, we have ω(0) = 0 and, using the Dominated Convergence Theorem (with majorizer sup s∈[0,1] ψ(s) < +∞ for small t ≥ 0), we conclude

ω + (0) = lim t 0 1 0 ϕ(st) st d s = lim t 0 1 0 ψ(st) d s = 1 0 lim t 0 ψ(st) d s = 0 .
It becomes a proper growth function, for example, if ψ(s) = 0 ⇔ s = 0 and we impose the additional condition lim t 0 w(t)/ψ(t) = 0. The function ψ(s) = cs α with α > 0, i.e. ϕ(s) = cs 1+α , is an example for a proper growth function.

Bregman distance. In order to introduce the notion of a Bregman function [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF], we first define a set of properties for functions to generate nicely behaving Bregman functions. Note that we have the duality (∂h) -1 = ∂h * where h * denotes the conjugate of h.

Definition 7 (Bregman distance [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF]). Let h : R N → R be proper, closed, convex and Gâteaux differentiable on int dom h = ∅. The Bregman distance associated with h is the function

D h : R N ×R N → [0, +∞] , (x, x) → h(x) -h(x) -x -x, ∇h(x) , if x ∈ int dom h ; +∞ , otherwise .
In contrast to the Euclidean distance, the Bregman distance is lacking symmetry.

We focus on Bregman distances that are generated by Legendre functions from the following class:

L :=      h : R N → R h is a proper, closed, convex Legendre function that is Fréchet differentiable on int dom h = ∅      .
To control the variable choice of Bregman distances throughout the algorithm's iterations, we introduce the following ordering relation for h 1 , h ∈ L :

h 1 h ⇔ ∀x ∈ dom h : ∀x ∈ int dom h : D h 1 (x, x) ≥ D h (x, x) .
As a consequence of h 1 h, we have dom

D h 1 ⊂ dom D h .
In order to conveniently work with Bregman distances, we collect a few properties.

Proposition 8. Let h ∈ L and D h be the associate Bregman distance.

(i) D h is strictly convex on every convex subset of dom ∂h with respect the first argument.

(ii) For x ∈ int dom h, it holds that D h (x, x) = 0 if and only if x = x.

(iii) For x ∈ R N and x, x ∈ int dom h the following three point identity holds:

D h (x, x) = D h (x, x) + D h (x, x) + x -x, ∇h(x) -∇h(x) .
Proof. (i) and (ii) follow directly from the definition of h being essentially strictly convex.

(iii) is stated in [START_REF] Chen | Convergence analysis of proximal-like minimization algorithm using bregman functions[END_REF]. It follows from the definition of a Bregman distance.

Associated with such a distance function is the following proximal mapping.

Definition 9 (Bregman proximal mapping [START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF]Def. 3.16]). Let f : R N → R and D h be a Bregman distance associated with h ∈ L . The D h -prox (or Bregman proximal mapping) associated with f is defined by

P h f (x) := argmin x f (x) + D h (x, x) . (1)
In general, the proximal mapping is set-valued, however for a convex function, the following lemma simplifies the situation. Lemma 10. Let f : R N → R be a proper, closed, convex function that is bounded from below, and h ∈ L such that int dom h ∩ dom f = ∅. Then the associated Bregman proximal mapping P h f is single-valued on its domain and maps to int dom h ∩ dom f . Proof. Single-valuedness follows from [START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF]Corollary 3.25(i)]. The second claim is from [START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF]Prop. 3.23(v)

(b)].
Proposition 11. Let f : R N → R be a proper, closed, convex function that is bounded from below, and h ∈ L such that int dom h ∩ dom f = ∅. For x ∈ int dom h, x = P h f (x), and any x ∈ dom f the following inequality holds:

f (x) + D h (x, x) ≥ f (x) + D h (x, x) + D h (x, x) .
Proof. See [START_REF] Chen | Convergence analysis of proximal-like minimization algorithm using bregman functions[END_REF]Lemma 3.2].

For examples and more useful properties of Bregman functions, we refer the reader to [START_REF] Bauschke | Legendre functions and the method of random Bregman projections[END_REF][START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF][START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF][START_REF] Nguyen | Forward-Backward Splitting with Bregman Distances[END_REF].

Miscellaneous. We make use of little-o notation f ∈ o(g) (or f = o(g)), which indicates that the asymptotic behavior of a function f is dominated by that of the function g. Formally, it is defined by

f ∈ o(g) ⇔ ∀ε > 0 : |f (x)| ≤ ε|g(x)| for |x| sufficiently small.
Note that a function ω is in o(t) if, and only if ω is a growth function.

Line Search Based Bregman Minimization Algorithms

In this paper, we solve optimization problems of the form min

x∈R N f (x) (2)
where f : R N → R is a proper, closed function on R N . We assume that Argmin f = ∅ and f := min f > -∞. The main goal is to develop a provably (subsequentially) convergent algorithm that finds a stationary point x of (2) in the sense of the limiting slope |∇f |(x) = 0. We analyze abstract algorithms that sequentially minimize convex models of the objective function.

The Abstract Algorithm

For each point x, we consider a proper, closed, convex model function

f x : R N → R with |f (x) -f x(x)| ≤ ω(|x -x|) , (3)
where ω is a growth function as defined in Definition 1. The model assumption (3) is an abstract description of a (local) first order oracle. For examples, we refer to Section 5.

Before delving further into the details, we need a bit of notation. Let

f h x,z (x) := f x(x) + D h (x, z) and f h x := f h x,x ,
where h ∈ L . Note that f h x (x) = f (x). Moreover, the following quantity defined for generic points x, x and x will be important:

∆ h x(x, x) := f h x (x) -f h x (x) . (4)
For x = x, it measures the decrease of the surrogate function f h x from the current iterate x to any other point x. Obviously, the definition implies that ∆ h

x(x, x) = 0 for all x.

Algorithm. We consider the following Algorithm 1.

Algorithm 1 (Inexact Bregman Proximal Minimization Line Search).

• Basic prerequisites: Fix γ ∈]0, 1[ and h ∈ L . Let

• (x k ) k∈N and (ỹ k ) k∈N be sequences in R N ;

• (f x k ) k∈N be a sequence of model functions with inf k∈N inf x f x k (x) > -∞;

• (h k ) k∈N be a sequence in L with h k h;

• (η k ) k∈N be a sequence of positive real numbers.

• Initialization: Select x 0 ∈ dom f ∩ int dom h.

• For each k ≥ 0: Generate the sequences such that the following relations hold:

∆ h k x k (ỹ k , x k ) < 0 with ỹk ∈ int dom h (5)
x k+1 = x k + η k (ỹ k -x k ) ∈ int dom h (6) f (x k+1 ) ≤ f (x k ) + γη k ∆ h k x k (ỹ k , x k ) (7)
If (5) cannot be satisfied, then the algorithm terminates.

The algorithm starts with a feasible point1 x 0 . At each iteration, it computes a point ỹk that satisfies [START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF], which is an inexact solution of the Bregman proximal mapping

xk = P h k fx k (x k ) := argmin x∈R N f x k (x) + D h k (x, x k ) (8)
that, at least, improves the (model) value compared to x k . Thanks to the class of Legendre functions L , this proximal mapping is well-defined and single-valued on its domain. The exact version of the algorithm solves the proximal mapping exactly for the global optimal solution. The optimal solution of the proximal mapping will always be denoted by xk instead of ỹk , which refers to an approximate solution. The direction ỹk -x k can be considered as a descent direction for the function f . Given this direction, the goal of ( 6) and ( 7) is the estimation of a step size η k (by line search, cf. Algorithm 2) that reduces the value of the objective function. In case that the proximal mapping has a solution but the first relation [START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF] can only be satisfied with equality, we will see that x k = ỹk must be a stationary point of the objective, hence, the algorithm terminates.

Remark 12. Instead of performing backtracking on the objective values as in Algorithm 1, backtracking on the scaling of the Bregman distance in ( 5) is also possible. For a special model function, this leads to ProxDescent [START_REF] Lewis | A proximal method for composite minimization[END_REF][START_REF] Drusvyatskiy | Error bounds, quadratic growth, and linear convergence of proximal methods[END_REF] (with Euclidean proximity function). If a scaled version of (5) yields a descent on f , we can set η k = 1, and accept this point. However, this can be expensive when the proximal subproblem in ( 5) is hard to solve, since each trial step requires to solve the subproblem. In order to break the backtracking, the new objective value must be computed anyway. Therefore, a computational advantage of the line search ( 6) and ( 7) is to be expected (cf. Section 6.1).

Algorithm 2 (Line Search for Algorithm 1).

• Basic prerequisites: Fix δ, γ ∈]0, 1[, η > 0, and k ∈ N.

• Input: Current iterates x k ∈ int dom h and ỹk satisfy (5).

• Solve: Find the smallest j ∈ N such that ηj := ηδ j satisfies ( 6) and ( 7).

• Return: Set the feasible step size η k for iteration k to ηj .

Algorithm 1-2 is well defined as the following lemmas show.

Lemma 13 (Well-definedness). Let ω in (3) be a growth function. Algorithm 1 is welldefined, i.e., for all k ∈ N, the following holds:

(i) there exists ỹk that satisfies (5) or x k = xk and the algorithm terminates;

(ii) x k ∈ dom f ∩ int dom h; and

(iii) there exists η k that satisfies ( 6) and [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF].

Proof. (i) For x k ∈ int dom h, Lemma 10 shows that P h k fx k maps to int dom h k ∩ dom f x k ⊂ int dom h ∩ dom f
and is single-valued. Thus, for example, ỹk = xk satisfies [START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF]. Otherwise, [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF] it holds that x k ∈ dom f for all k. Since x k ∈ int dom h and ỹk ∈ dom h, for small η k also x k+1 ∈ int dom h, hence x k+1 ∈ dom f ∩ int dom h. Inductively, we conclude the statement. (iii) This will be shown in Lemma 14.

x k = xk , which shows (i). (ii) Since x 0 ∈ dom f ∩ int dom h and f (x k+1 ) ≤ f (x k ) by
Lemma 14 (Finite termination of Algorithm 2). Consider Algorithm 1 and fix k ∈ N. Let ω in (3) be a growth function. Let δ, γ ∈]0, 1[, η > 0, h := h k , and x := x k , ỹ := ỹk be such that ∆ h x(ỹ, x) < 0. Then, there exists j ∈ N such that ηj := ηδ j satisfies

f (x + ηj (ỹ -x)) ≤ f (x) + γ ηj ∆ h x(ỹ, x) .
Proof. This result is proved by contradiction. Define v := ỹ -x. By our assumption in (3), we observe that

f (x + ηj v) -f (x) ≤ f x(x + ηj v) -f (x) + o(η j ) . (9)
Using Jensen's inequality for the convex function f x provides:

f x(x + ηj v) -f x(x) ≤ ηj f x(x + v) + (1 -ηj )f x(x) -f x(x) = ηj • f x(x + v) -f x(x) . (10)
Now, suppose γ∆ h x(ỹ, x) < 1 ηj (f (x + ηj v) -f (x)) holds for any j ∈ N. Then, using ( 9) and ( 10), we conclude the following:

γ∆ h x(ỹ, x) < f x(x + v) -f x(x) + o(η j )/η j ≤ f x(x + v) -f x(x) + Dh(ỹ, x) + o(η j )/η j = (f h x (ỹ) -f h x (x)) + o(η j )/η j = ∆ h x(ỹ, x) + o(η j )/η j ,
which for j → ∞ yields the desired contradiction, since γ ∈]0, 1[ and ∆ h x(ỹ, x) < 0.

Finite Time Convergence Analysis

First, we study the case when the algorithm terminates after a finite number of iterations, i.e., there exists k 0 ∈ N such that (5) cannot be satisfied. Then, the point ỹk 0 is a global minimizer of f

h k 0 x k 0 and ∆ h k 0 x k 0 (ỹ k 0 , x k 0 ) = 0.
Moreover, the point x k 0 turns out to be a stationary point of f . Lemma 15. For x ∈ dom f and a model f x that satisfies (3), where ω is a growth function, the following holds:

|∇f x|(x) = |∇f |(x) .
Proof. Since ω(0) = 0, we have from (3) that f x(x) = f (x). This, together with sub-

additivity of [•] + , entails [f x(x) -f x(x)] + |x -x| ≤ [f (x) -f (x)] + + [f (x) -f x(x)] + |x -x| ≤ [f (x) -f (x)] + |x -x| + |f (x) -f x(x)| |x -x| ≤ [f (x) -f (x)] + |x -x| + ω(|x -x|) |x -x|
Passing to the lim sup on both sides and using that ω ∈ o(t), we get

|∇f x|(x) ≤ |∇f |(x).
Arguing similarly but now starting with |∇f |(x), we get the reverse inequality, which in turn shows the claimed equality.

Proposition 16 (Stationarity for finite time termination). Consider the setting of Algorithm 1. Let ω in (3) be a growth function. Let k 0 ∈ N be fixed, and set x = ỹk 0 , x = x k 0 , h = h k 0 , and x, x ∈ dom f ∩ int dom h. If ∆ h x(x, x) ≥ 0, then x = x, ∆ h x(x, x) = 0, and |∇f |(x) = 0, i.e. x is a stationary point of f .

Proof. Since x is the unique solution of the proximal mapping, obviously ∆ h x(x, x) = 0 and x = x. Moreover, x is the minimizer of f h x , i.e. we have

0 = |∇f h x |(x) = |∇f h x |(x) = lim sup x→x x =x [f x(x) -f x(x) -Dh(x, x)] + |x -x| = |∇f x|(x) = |∇f |(x) ,
where we used that h is Fréchet differentiable at x and Lemma 15.

Asymptotic Convergence Analysis

We have established stationarity of the algorithm's output, when it terminates after a finite number of iterations. Therefore, without loss of generality, we now focus on the case where ( 5) can be satisfied for all k ∈ N. We need to make the following assumptions.

Assumption 1. The sequence (ỹ k ) k∈N satisfies f h k x k (ỹ k ) ≤ inf f h k x k + ε k for some ε k → 0.
Remark 17. Assumption 1 states that asymptotically (for k → ∞) the Bregman proximal mapping ( 8) must be solved accurately. In order to obtain stationarity of a limit point, Assumption 1 is necessary, as shown by Bonettini et al. [START_REF] Bonettini | Variable metric inexact line-search based methods for nonsmooth optimization[END_REF]after Theorem 4.1] for a special setting of model functions.

Assumption 2. Let h ∈ L . For every bounded sequences (x k ) k∈N and (x k ) k∈N in int dom h, and (h k ) k∈N such that h k h, it is assumed that:

x k -xk → 0 ⇔ D h k (x k , xk ) → 0 .
Remark 18. (i) Assumption 2 states that (asymptotically) a vanishing Bregman distance reflects a vanishing Euclidean distance. This is a natural assumption and satisfied, e.g., by most entropies such as Boltzmann-Shannon, Fermi-Dirac, and Burg entropy.

(ii) The equivalence in Assumption 2 is satisfied, for example, when there exists c ∈ R such that c h h k holds for all k ∈ N and the following holds:

x k -xk → 0 ⇔ D h (x k , xk ) → 0 .
Proposition 19 (Convergence of objective values). Consider the setting of Algorithm 1.

Let ω in (3) be a growth function. The sequence of objective values (f (x k )) k∈N is nonincreasing and converging to some

f * ≥ f > -∞.
Proof. This statement is a consequence of ( 7) and ( 5), and the lower-boundedness of f .

Asymptotically, under some condition on the step size, the improvement of the model objective value between ỹk and x k must tend to zero. Since we do not assume that the step sizes η k are bounded away from zero, this is a non-trivial result.

Proposition 20 (Vanishing model improvement). Consider the setting of Algorithm 1. Let ω in (3) be a growth function. Suppose, either2 inf k η k > 0 or η k is selected by the Line Search Algorithm 2. Then,

∞ k=0 η k (-∆ h k x k (ỹ k , x k )) < +∞ and ∆ h k x k (ỹ k , x k ) → 0 as k → ∞.
Proof. The first part follows from rearranging [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF], and summing both sides for k = 0, . . . , n:

γ n k=0 η k (-∆ h k x k (ỹ k , x k )) ≤ n k=0 (f (x k ) -f (x k+1 )) = f (x 0 ) -f (x n+1 ) ≤ f (x 0 ) -f * .
In the remainder of the proof, we show that ∆ h k x k (ỹ k , x k ) → 0, which is not obvious unless inf k η k > 0. The model improvement is bounded. Boundedness from above is satisfied by construction of the sequence (ỹ k ) k∈N . Boundedness from below follows from the following observation and the uniform boundedness of the model functions from below:

∆ h k x k (ỹ k , x k ) = f h k x k (ỹ k ) -f h k x k (x k ) ≥ f h k x k (x k ) -f (x k ) ≥ f x k (x k ) -f (x 0 ) . Therefore, there exists K ⊂ N such that the subsequence ∆ h k x k (ỹ k , x k ) converges to some ∆ * as k K → ∞. Suppose ∆ * < 0.
Then, the first part of the statement implies that the step size sequence must tend to zero, i.e., η k → 0 for k K → ∞. For k ∈ K sufficiently large, the line search procedure in Algorithm 2 reduces the step length from η k /δ to η k . (Note that η k can be assumed to be the "first" step length that achieves a reduction in [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF]). Before multiplying with δ, no descent of (7) was observed, i.e.,

(η k /δ)γ∆ h k x k (ỹ k , x k ) < f (x k + (η k /δ)v k ) -f (x k )
, where v k = ỹk -x k . Using ( 9) and [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF], we can make the same observation as in the proof of Lemma 14:

γ∆ h k x k (ỹ k , x k ) < f x k (x k + v) -f x k (x k ) + o(η k /δ)/(η k /δ) ≤ f x k (x k + v) -f x k (x k ) + D h k (ỹ k , x k ) + o(η k )/η k = (f h k x k (ỹ k ) -f h k x k (x k )) + o(η k )/η k = ∆ h k x k (ỹ k , x k ) + o(η k )/η k , which for η k → 0 yields a contradiction, since γ ∈]0, 1[ and ∆ h k x k (ỹ k , x k ) < 0. Therefore, any cluster point ∆ * of (∆ h k x k (ỹ k , x k
)) k∈K must be 0, which concludes the proof.

Asymptotic Stationarity with a Growth Function

In order to establish stationarity of limit points generated by Algorithm 1 additional assumptions are required. We consider three different settings for the model assumption [START_REF] Bauschke | Legendre functions and the method of random Bregman projections[END_REF]: ω in the model assumption (3) is a growth function (this section), ω is a proper growth function (Section 4.3.2), and ω is global growth function of the form ω = D h (Section 4.3.3).

Assumption 3. Let x * be a limit point of (x k ) k∈N and Using this assumption, we can state one of our main theorems, which shows convergence to a stationary point under various condition. The conditions are easily verified in many applications (see Section 5).

x k f → x * as k K → ∞ with K ⊂ N. Then |∇f x k |(x k ) = |∇f |(x k ) → 0 as k K → ∞ .
Theorem 22 (Asymptotic stationarity with a growth function). Consider the setting of Algorithm 1. Let ω in (3) be a growth function. Moreover, let either inf k η k > 0 or η k be selected by the Line Search Algorithm 2. Let (x k ) k∈N and (ỹ k ) k∈N be bounded sequences such that Assumptions 1 and 2 hold and let f x k obey (3) with growth function ω. Then,

x k -ỹk → 0 and for xk = P h k fx k (x k ), it holds that x k -xk → 0 and xk -ỹk → 0. Moreover, f (x k ) -f (ỹ k ) → 0 and f (x k ) -f (x k ) → 0 as k → ∞. Suppose Assumption 3 is satisfied. If x *
is a limit point of the sequence (x k ) k∈N , and one of the following conditions is satisfied:

(i) f is continuous on the closure of dom h, (ii) x * ∈ int dom h, (iii) x * ∈ dom h and D h k (x * , ỹk ) → 0 as k K → ∞, (iv) x * ∈ cl dom h and • for all x ∈ int dom h ∩ dom f holds that D h k (x, xk ) -D h k (x, x k ) → 0 as k K → ∞,
• and for all x ∈ dom f the model functions obey

f x k (x) → f x * (x) as k K → ∞, then x * is a stationary point of f .
Proof. First, we show that for k → ∞ the pairwise distances between the sequences (x k ) k∈N , (ỹ k ) k∈N , and (x k ) k∈N vanishes. Proposition 11, reformulated in our notation, can be stated as

∆ h k x k (x, xk ) = f h k x k (x) -f h k x k (x k ) ≥ D h k (x, xk ) , ∀x ∈ dom f . (11)
As a direct consequence, using x = ỹk together with Assumptions 1 and 2, we obtain

D h k (ỹ k , xk ) → 0 thus xk -ỹk → 0 .
Moreover, from Proposition 20, we have ∆ h k x k (ỹ k , x k ) → 0, and from

∆ h k x k (ỹ k , x k ) = ∆ h k x k (ỹ k , xk ) -∆ h k x k (x k , xk ) ≤ ∆ h k x k (ỹ k , xk ) -D h k (x k , xk ) , (12)
and Assumptions 1 and 2, we conclude that x k -xk → 0, hence also x k -ỹk → 0.

The next step is to show that f (x k ) -f (ỹ k ) → 0 as k → ∞. This follows from the following estimation:

|f (x k ) -f (ỹ k )| ≤ |f x k (x k ) -f x k (ỹ k )| + ω(|ỹ k -x k |) ≤ |f h k x k (x k ) -f h k x k (ỹ k )| + D h k (ỹ k , x k ) + ω(|ỹ k -x k |) = |∆ h k x k (x k , ỹk )| + D h k (ỹ k , x k ) + ω(|ỹ k -x k |) , (13)
where the right hand side vanishes for k → ∞. Analogously, we can show that f

(x k ) - f (x k ) → 0 as k → ∞.
Let x * be the limit point of the subsequence (x k ) k∈K for some K ⊂ N. The remainder of the proof shows that f (i) implies f (ỹ k ) → f (x * ) as k → ∞ by definition. For (ii) and (iii), we make the following observation:

(ỹ k ) → f (x * ) as k → ∞. Then f (x k ) -f (ỹ k ) → 0 implies that x k f → x * as k K → ∞,
f (ỹ k )-ω(|ỹ k -x k |) ≤ f h k x k (ỹ k ) = f h k x k (x k )+(f h k x k (ỹ k )-f h k x k (x k )) ≤ f h k x k (x * )+∆ h k x k (ỹ k , xk ) , (14)
where xk = P h k fx k (x k ). Taking "lim sup k K →∞ " on both sides, D h k (x * , x k ) → 0 (Assumption 2 for (ii) or the assumption in (iii)), and ∆ h k

x k (ỹ k , xk ) → 0 (Assumption 1) shows that lim sup

k K →∞ f (ỹ k ) ≤ f (x * ). Since f is closed, f (ỹ k ) → f (x * ) holds.
We consider (iv). For all x ∈ int dom h ∩ dom f , we have [START_REF] Bonettini | Variable metric inexact line-search based methods for nonsmooth optimization[END_REF] or, reformulated,

f h k x k (x) - D h k (x, xk ) ≥ f h k
x k (x k ), which implies the following:

f x k (x) + D h k (x, x k ) -D h k (x, xk ) -D h k (x k , x k ) ≥ f (x k ) -ω(|x k -x k |) .
Note that for any x the limits for k K → ∞ on the left hand side exist. In particular, we have

D h k (x, x k ) -D h k (x, xk ) -D h k (x k , x k ) → 0 as k K → ∞ ,
by the assumption in (iv), and Assumption 2 together with xk -x k → 0. The limit of f x k (x) exists by assumption and coincides with

f x * (x). Choosing a sequence (z k ) k∈N in int dom h ∩ dom f with z k → x * as k K → ∞, in the limit, we obtain f (x * ) ≥ lim k K →∞ f (x k ) =: f * , since f x * (z k ) → f x * (x * ) = f (x * ) for z k → x * as k K → ∞. Invoking that f is closed, we conclude the f -attentive convergence f * = lim inf k→∞ f (x k ) ≥ f (x * ) ≥ f * .
Remark 23. Existence of a limit point x * is guaranteed by assuming that (x k ) k∈N is bounded. Alternatively, we could require that f is coercive (i.e. f (x k ) → ∞ for |x k | → ∞), which implies boundedness of the lower level sets of f , hence by Proposition 19 the boundedness of (x k ) k∈N .

Remark 24. From Theorem 22, clearly, also ỹk

f → x * and xk f → x * as k K → ∞ holds.
Therefore, Assumption 3 could also be stated as the requirement

|∇f |(x k ) → 0 or |∇f |(ỹ k ) → 0 as k K → ∞ ,
in order to conclude that limit points of (x k ) k∈N are stationary points.

As a simple corollary of this theorem, we replace Assumption 3 with the relative error condition mentioned in Remark 21.

Corollary 25 (Asymptotic stationarity with a growth function). Consider the setting of Algorithm 1. Let ω in (3) be a growth function. Moreover, let either inf k η k > 0 or η k be selected by the Line Search Algorithm 2. Let (x k ) k∈N and (ỹ k ) k∈N be bounded sequences such that Assumptions 1 and 2 hold and let f x k obey (3) with growth function ω. Suppose there exists a continuous function ψ :

R + → R + with ψ(0) = 0 such that |∇f |(x k+1 ) ≤ ψ(|x k+1 -x k |) is satisfied. If x *
is a limit point of the sequence (x k ) k∈N and one of the conditions (i)-(iv) in Theorem 22 is satisfied, then x * is a stationary point of f . Proof. Theorem 22 shows that ỹk -x k → 0, thus, inf k η k > 0 implies x k+1 -x k → 0 by [START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF]. Therefore, the relation |∇f |(x k+1 ) ≤ ψ(|x k+1 -x k |) shows that Assumption 3 is automatically satisfied and we can apply Theorem 22 to deduce the statement. Some more results on the limit point set. In Theorem 22 we have shown that limit points of the sequence (x k ) k∈N generated by Algorithm 1 are stationary, and in fact the sequence f -converges to its limit points. The following proposition shows some more properties of the set of limit points of (x k ) k∈N . This is a well-known result [10, Lem. 5] that follows from x k+1 -x k → 0 as k → ∞. Proposition 26. Consider the setting of Algorithm 1. Let ω in (3) be a growth function and inf k η k > 0. Let (x k ) k∈N and (ỹ k ) k∈N be bounded sequences such that Assumptions 1, 2 and 3 hold. Suppose one of the conditions (i)-(iv) in Theorem 22 is satisfied for each limit point of (x k ) k∈N . Then, the set

S := x * ∈ R N | ∃K ⊂ N : x k → x * as k K → ∞ of limit points of (x k ) k∈N is connected, each point x * ∈ S is stationary for f , and f is constant on S.
Proof. Theorem 22 shows that ỹk -x k → 0. Thus, boundedness of η k away from 0 implies x k+1 -x k → 0 by [START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF]. Now, the statement follows from [10, Lem. 5] and Theorem 22.

Asymptotic Stationarity with a Proper Growth Function

Our proof of stationarity of limit points generated by Algorithm 1 under the assumption of a proper growth function ω in (3) relies on an adaptation of a recently proved result by Drusvyatskiy et al. [START_REF] Drusvyatskiy | Nonsmooth optimization using Taylor-like models: error bounds, convergence, and termination criteria[END_REF]Corollary 5.3], which is stated in Lemma 27 before the main theorem of this subsection. The credits for this lemma should go to [START_REF] Drusvyatskiy | Nonsmooth optimization using Taylor-like models: error bounds, convergence, and termination criteria[END_REF].

Lemma 27 (Perturbation result under approximate optimality). Let f : R N → R be a proper closed function. Consider bounded sequences (x k ) k∈N and (ỹ k ) k∈N with x k -ỹk → 0 for k → ∞, and model functions f h k

x k according to (3) with proper growth functions. Suppose Assumption 1 and 2 hold. If (x * , f (x * )) is a limit point of (x k , f (x k )) k∈N , then x * is stationary for f . Proof. Recall ε k from Assumption 1. Theorem 5.1 from [START_REF] Drusvyatskiy | Nonsmooth optimization using Taylor-like models: error bounds, convergence, and termination criteria[END_REF] guarantees, for each k and any ρ k > 0, the existence of points ŷk and z k such that the following hold:

(i) (point proximity) |ỹ k -z k | ≤ ε k ρ k and |z k -ŷk | ≤ 2 • ω(|z k -x k |) ω + (|z k -x k |) , under the convention 0 0 = 0 , (ii) (value proximity) f (ŷ k ) ≤ f (ỹ k ) + 2ω(|z k -x k |) + ω(|ỹ k -x k |), and (iii) (near-stationarity) |∇f |(ŷ) ≤ ρ k + ω + (|z k -x k |) + ω + (|ŷ k -x k |), Setting ρ k = √ ε k , using ε k → 0 and the point proximity, shows that |ỹ k -z k | → 0. Moreover |z k -x k | ≤ |z k -ỹk | + |ỹ k -x k | → 0, which implies that |z k -ŷk | → 0. Now, we fix a convergent subsequence (x k , f (x k )) → (x * , f (x * )) as k K → ∞ for some K ⊂ N.
Using (11), we observe xk -ỹk → 0, hence x k -xk → 0. From Proposition 20 and Assumption 1, we conclude that ∆ h k x k (ỹ k , x k ) → 0, and, therefore f (x k ) -f (ỹ k ) → 0 using [START_REF] Burg | The relationship between maximum entropy spectra and maximum likelihood spectra[END_REF]. Consider the value proximity. Combined with the lower semi-continuity of f , it yields

f (x * ) ≤ lim inf k K →∞ f (ŷ k ) ≤ lim sup k K →∞ f (ŷ k ) ≤ lim sup k K →∞ f (ỹ k ) ≤ f (x * ) , hence (ŷ k , f (ŷ k )) → (x * , f (x * )) as k K → ∞. Near-stationarity implies that |∇f |(ŷ k ) → 0, which proves that |∇f |(x * ) = 0, hence x * is a stationary point.
Remark 28. The setting in [START_REF] Drusvyatskiy | Nonsmooth optimization using Taylor-like models: error bounds, convergence, and termination criteria[END_REF]Corollary 5.3] is recovered when (x k ) k∈N is given by x k+1 = ỹk .

Theorem 29 (Asymptotic stationarity with a proper growth function). Consider the setting of Algorithm 1. Let ω in (3) be a proper growth function. Moreover, let either inf k η k > 0 or η k be selected by the Line Search Algorithm 2. Let (x k ) k∈N and (ỹ k ) k∈N be bounded sequences such that Assumptions 1 and 2 hold. If x * is a limit point of the sequence x k and one of the conditions (i)-(iv) in Theorem 22 is satisfied, then x * is a stationary point of f . Proof. Propositions 19 and 20, and the proof of f -attentive convergence from Theorem 22 only rely on a growth function. Instead of assuming that the slope vanishes, here we apply Lemma 27 to conclude stationarity of the limit points.

Of course, Proposition 26 can also be stated in the context here.

Asymptotic Analysis with a Global Growth Function

Suppose, for x ∈ int dom h for some h ∈ L , the model error can be estimated as follows:

|f (x) -f x(x)| ≤ LD h (x, x) ∀x . ( 15
)
Since h is Fréchet differentiable on int dom h, the right hand side is bounded by a growth function. Without loss of generality, we restrict ourselves to a fixed function h ∈ L (this section analyses a single iteration). In order to reveal similarities to well-known step size rules, we scale h in the definition of f h x to D h/α = 1 α D h with α > 0 instead of D h . Here, decreasing objective values can be assured without the line search procedure (see Proposition 30), i.e., η k = 1 is always feasible.

In order to obtain the result of stationarity of limit points (Theorem 22 or 29), we can either verify by hand that Assumption 3 holds or we need to assume that D h (x, x) is bounded by a proper growth function.

Proposition 30. Consider the setting of Algorithm 1 and let (15) be satisfied.

(i) For points ỹ that satisfy ∆

h x(ỹ, x) < 0, 1 -αL α D h (ỹ, x) ≤ f (x) -f (ỹ)
holds, where the left-hand-side is strictly larger than 0 for α ∈]0, 1/L[.

(ii) For points x = P h fx (x), the following descent property holds:

1 + ρ -αL α D h (x, x) ≤ f (x) -f (x) ,
where the left-hand-side is strictly larger than 0 for α ∈]0, (1 + ρ)/L[, and ρ is the Bregman symmetry factor defined by ρ := inf D h (x,x) D h (x,x) | x, x ∈ int dom h , x = x ; (see [START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF]).

Proof. The following relations hold:

∆ h x(ỹ, x) ≤ 0 ⇔ f h x (ỹ) ≤ f h x (x) ⇔ f x(ỹ) + 1 α D h (ỹ, x) ≤ f x(x) = f (x) . (16)
Bounding the left hand side of the last expression using [START_REF] Chaudhuri | Blind Image Deconvolution[END_REF], we obtain

f (ỹ) -LD h (ỹ, x) + 1 α D h (ỹ, x) ≤ f (x) , (17)
which proves part (i). Part (ii) follows analogously. However, thanks to the three point inequality from Proposition 11 and optimality of x the rightmost inequality of ( 16) improves to

f x(x) + 1 α D h (x, x) + 1 α D h (x, x) ≤ f x(x) = f (x) ,
and the statement follows.

A Remark on Convex Optimization

In this section, let f be convex, and consider the following global model assumption

0 ≤ f (x) -f x(x) ≤ LD h (x, x) . ( 18
)
We establish a convergence rate of O(1/k) for Algorithm 1 with η k ≡ 1. For Forward-Backward Splitting, this has been shown by Bauschke et al. [START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF]. We only require f x to be a model w.r.t. [START_REF] Combettes | Dualization of signal recovery problems[END_REF].

Proposition 31. Consider Algorithm 1 with η k ≡ 1 and model functions that obey [START_REF] Combettes | Dualization of signal recovery problems[END_REF]. For x k+1 = P h/α fx k (x k ) and α = 1 L , the following rate of convergence on the objective values holds:

f (x k+1 ) -f (x) ≤ LD h (x * , x 0 ) 2k (= O(1/k)) .
Proof. The three point inequality in Proposition 11 combined with the model assumption (18) yields the following inequality:

f (x) + 1 -αL α D h (x, x) + 1 α D h (x, x) ≤ f (x) + 1 α D h (x, x)
for all x. Restricting to 0 < α ≤ 1 L , we obtain

f (x) -f (x) ≤ 1 α (D h (x, x) -D h (x, x)) . (19)
Let x * be a minimizer of f . We make the following choices:

x = x * , x = x k+1 , and x = x k .

Summing both sides up to iteration k and the descent property yield the convergence rate:

f (x k+1 ) -f (x) ≤ D h (x * , x 0 ) 2αk α= 1 L = LD h (x * , x 0 ) 2k . ( 20
)
Example 32. This section goes beyond the Forward-Backward Splitting setting, e.g., we may specialize Example 35 to convex problems of the form min x g(F (x)), for instance, using

F (x) = (f 1 (x), f 2 (x)) and g(z 1 , z 2 ) = max{z 1 , z 2 }.
We discuss several classes of problems that can be solved using our framework. To apply Algorithm 1, in Section 5.1, we define a suitable model and mention the associated algorithmic step that arises from exactly minimizing the sum of the model and an Euclidean proximity measure. However our algorithm allows for inexact solutions and very flexible (also iteration dependent) Bregman proximity functions. Examples are provided in Section 5.2. For brevity, we define the symbols Γ 0 for the set of proper, closed, convex functions and C 1 for the set of continuously differentiable functions.

Examples of Model Functions

Example 33 (Forward-Backward Splitting). Problems of the form

f = f 0 + f 1 with f 0 ∈ Γ 0 and f 1 ∈ C 1 can be modeled by f x(x) = f 0 (x) + f 1 (x) + x -x, ∇f 1 (x) .
This model is associated with Forward-Backward Splitting (FBS). We assume that one of the following error models is satisfied:

|f (x)-f x(x)| = |f 1 (x)-f 1 (x)-x -x, ∇f 1 (x) | ≤          L 2 |x -x| 2 , if ∇f 1 is L-Lipschitz ; 1 0 ϕ(t|x-x|) t d t , if ∇f 1 is ψ-uniformly continuous ; ω(|x -x|) , otherwise ,
which is the linearization error of the smooth part f 1 . The first case obeys a global (proper) growth function, derived from the common Descent Lemma. The second case is the generalization to a ψ-uniformly continuous gradient with ϕ(s) = sψ(s) as in Lemma 4. The bound

1 0 ϕ(t|x-x|) t
d t is a growth function but not necessarily a proper growth function. The third case is the most general and assumes that the error obeys a growth function. In any case, the model satisfies the model consistency required in Theorem 22(iv). For any x ∈ dom f and x → x * ,

|f 0 (x) + f 1 (x) + x -x, ∇f 1 (x) -(f 0 (x) + f 1 (x * ) + x -x * , ∇f 1 (x * ) )| → 0
holds, thanks to the continuous differentiability of f 1 and continuity of the inner product.

In order to verify Assumption 3, we make use of Remark 24 and show that

|∇f |(x k ) → 0 as k K → ∞ where K ⊂ N is such that x k K → x * .
Note that xk satisfies the following relation:

0 ∈ ∂f 0 (x k ) + ∇f 1 (x k ) + ∇h k (x k ) -∇h k (x k ) ⇒ ∇f 1 (x k ) -∇f 1 (x k ) + ∇h k (x k ) -∇h k (x k ) ∈ ∂f 0 (x k ) + ∇f 1 (x k ) = ∂f (x k ) Moreover, we know that xk -x k → 0 as k K → ∞. Since ∇f 1 is continuous, if |∇h k (x k ) - ∇h k (x k )| → 0 for k K → ∞, then Assumption 3/Remark 24 is satisfied. The condition |∇h k (x k ) -∇h k (x k )| → 0 is
naturally fulfilled by many Legendre functions, e.g., if ∇h k is ψ-uniformly continuous (uniformly in k) with α > 0 or uniformly continuous (independent of k) on bounded sets or continuous at x * (uniformly w.r.t. k), and will be discussed in more detail in Section 5.2.

Example 34 (Variable metric FBS). We consider an extension of Examples 33. An alternative feasible model for a twice continuously differentiable function f 1 is the following:

f x(x) = f 0 (x) + f 1 (x) + x -x, ∇f 1 (x) + 1 2 x -x, B(x -x) ,
where B := [∇ 2 f 1 (x)] + is a positive definite approximation to ∇ 2 f 1 (x), which leads to a Hessian driven variable metric FBS. It is easy to see that the model error satisfies the growth function ω(s). Again, Theorem 22(iv) obviously holds and the same conclusions about Assumption 3 can be made as in Example 33.

Example 35 (ProxDescent). Problems of the form

f 0 + g • F with f 0 ∈ Γ 0 , F ∈ C 1 , and g ∈ Γ 0 finite-valued ,
which often arise from non-linear inverse problems, can be approached by the model function

f x(x) = f 0 (x) + g(F (x) + DF (x)(x -x)) ,
where DF (x) is the Jacobian matrix of F at x. The associated algorithm is connected to ProxDescent [START_REF] Lewis | A proximal method for composite minimization[END_REF][START_REF] Drusvyatskiy | Error bounds, quadratic growth, and linear convergence of proximal methods[END_REF]. If g is a quadratic function, the algorithm reduces to the Levenberg-Marquardt algorithm [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF]. The error model can be computed as follows:

|f (x) -f x(x)| = |g(F (x)) -g(F (x) + DF (x)(x -x))| ≤ |F (x) -F (x) -DF (x)(x -x)| ≤        L 2 |x -x| 2 , if DF is L-Lipschitz and g is -Lipschitz ; 1 0 ϕ(t|x-x|) t d t , if DF is ψ-uniform continuous and g is -Lipschitz ; ω(|x -x|) , otherwise , (21)
where is the (possibly local) Lipschitz constant of g around F (x). Since g is convex and finite-valued, it is always locally Lipschitz continuous. Since F is continuously differentiable, for x sufficiently close to x, both F (x) and F (x)+DF (x)(x-x) lie in a neighborhood of F (x) where the local Lipschitz constant of g is valid, which shows the first inequality in [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF].

The second case uses the concept of ψ-uniform continuity from Definition 2 and Lemma 4. The third case is the assumption that the error obeys a growth function. With a similar reasoning, we can show that Theorem 22(iv) is satisfied.

We consider Assumption 3 (see also Remark 24). Let x k → x * as k K → ∞ for K ⊂ N and xk -x k → 0. Since g is finite-valued, using [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF]Corollary 16.38] (sum-rule for the subdifferential), and [START_REF] Rockafellar | Variational Analysis[END_REF]Theorem 10.6], we observe that

0 ∈ ∂f 0 (x k ) + DF (x k ) * ∂g(F (x k ) + DF (x k )(x k -x k )) + ∇h k (x k ) -∇h k (x k ) , (22)
where DF (x k ) * denotes the adjoint of DF (x k ). We can assume that, for k large enough, F (x k )+DF (x k )(x k -x k ) and F (x k ) lie a neighborhood of F (x * ) on which g has the Lipschitz constant > 0. By [START_REF] Rockafellar | Variational Analysis[END_REF]Theorem 9.13], ∂g is locally bounded around F (x * ), i.e. there exists a compact set G such that ∂g(z) ⊂ G for all z in a neighborhood of F (x * ). We conclude that sup

v∈∂g(F (x k )+DF (x k )(x k -x k )) w∈∂g(F (x k )) |DF (x k ) * v -DF (x k ) * w| ≤ sup v,w∈G |DF (x k ) * v -DF (x k ) * w| → 0 for k K → ∞ since DF (x k ) → DF (x * ) and DF (x k ) → DF (x *
). Again assuming that ∇h k (x k ) -∇h k (x k ) → 0 we conclude that the outer set-limit of the right hand side of ( 22) is included in ∂f (x k ) and, therefore, the slope |∇f |(x k ) vanishes for k

K → ∞.
Example 36. Problems of the form

f 0 +g•F with f 0 ∈ Γ 0 , g ∈ C 1
, ∇g i ≥ 0 , and F = (F 1 , . . . , F M ) is Lipschitz with F i ∈ Γ 0 can be modeled by

f x(x) = f 0 (x) + g(F (x)) + F (x) -F (x), ∇g(F (x)) .
Such problems appear for example in non-convex regularized imaging problems in the context of iteratively reweighted algorithms [START_REF] Ochs | On iteratively reweighted algorithms for nonsmooth nonconvex optimization in computer vision[END_REF]. For the error of this model function, we observe the following:

|f (x) -f x(x)| = |g(F (x)) -(g(F (x)) + F (x) -F (x), ∇g(F (x)) )| ≤      2 |F (x) -F (x)| 2 , if ∇g is -Lipschitz ; 1 0 ϕ(t|F (x)-F (x)|) t d t , if ∇g is ψ-uniform continuous ; ω(|F (x) -F (x)|) , otherwise ; ≤          L 2 2 |x -x| 2 , if ∇g is -Lipschitz and F is L-Lipschitz ; 1 0 ϕ(tL|x-x|)) t d t , if ∇g is ψ-uniformly continuous and F is L-Lipschitz ; ω(|x -x|) , otherwise ,
which shows the same growth functions are obeyed as in Example 33 and 35. The explanation for the validity of the reformulations are analogue to those of Example 35. It is easy to see that Theorem 22(iv) holds.

We consider Assumption 3/Remark 24. Let x k → x * as k K → ∞ for K ⊂ N and xkx k → 0. Since g is continuously differentiable, the sum-formula for the subdifferential holds. Moreover, we can apply [START_REF] Rockafellar | Variational Analysis[END_REF]Corollary 10.09] (addition of functions) to see that xk satisfies the following relation:

0 ∈ ∂f 0 (x k ) + M i=1 ∂F i (x k )(∇g(F (x k ))) i + ∇h k (x k ) -∇h k (x k ) ,
Note that by [START_REF] Rockafellar | Variational Analysis[END_REF]Theorem 10.49] the subdifferential of g

•F at xk is M i=1 ∂F i (x k )(∇g(F (x k ))) i .
As in Example 35, using the Lipschitz continuity of F , hence local boundedness of ∂F , and using the continuous differentiability of g, the sequence of sets

M i=1 ∂F i (x k )(∇g(F (x k ))) i - ∂F i (x k )(∇g(F (x k ))) i vanishes for k K → ∞, which implies that the slope |∇f |(x k ) vanishes for k K → ∞.
Example 37 (Problem adaptive model function). Our framework allows for a problem specific adaptation using a combination of Examples 33,[START_REF] Olshausen | Sparse coding with an overcomplete basis set: A strategy employed by V1?[END_REF][START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF][START_REF] Rockafellar | Variational Analysis[END_REF]. Consider the following objective function 3

f : R N × R 2 → R with a ∈ R N and b ∈ R: f (x, z) = f 1 (z) + δ [-1,1] 2 (z) + f 2 (x) with f 1 ∈ C 2 (R 2 ) is strongly convex ; f 2 (x) := max{( a, x -b) 2 , 1 -exp(-|x|)} .
We define our model function as:

f (x,z) (x, z) = f1 (z; z) + δ [-1,1] (z) + f2 (x; x) with f1 (z; z) := f 1 (z) + ∇f 1 (z), z -z + 1 2 z -z, ∇ 2 f 1 (z)(z -z) ; f2 (x; x) := max{( a i , x -b i ) 2 , 1 + exp(-|x|)(|x| -|x| -1)} .
The strong convexity of f 1 allows for a convex second order approximation with positive definite Hessian. We linearize only the second component of the "max" (w.r.t. |x|) in f 2 to obtain a convex approximation that is as close as possible to original function f 2 .

As ∇f 1 is Lipschitz continuous on the compact set [0, 1] 2 , the growth function w.r.t. z is of the form L|z -z| 2 . Moreover, using 1-Lipschitz continuity of exp on R -, we can bound 2 and, using Lipschitz continuity of |x|, by |x -x| 2 . Therefore, the model error is given by a growth function w(t) = max{1, L}t 2 .

| exp(-|x|) -exp(-|x|)(1 + |x| -|x|)| by ||x| -|x||

Examples of Bregman functions

Let us explore some of the Bregman functions, that are most important to our applications and show that our assumptions are satisfied.

Example 38 (Euclidean Distance). The most natural Bregman proximity function is the Euclidean distance

D h (x, x) = 1 2 |x -x| 2 ,
which is generated by the Legendre function h(x) = 1 2 |x| 2 . The domain of h is the whole space R N , which implies that Condition (ii) in Theorem 22 is satisfied for any limit point. Assumption 2 is trivial, and for the model functions in Section 5.1, Assumption 3 is satisfied, if x k -xk → 0 implies ∇h(x k ) -∇h(x k ) → 0, which is clearly true. Therefore, for the models in Section 5.1 combined with the Euclidean proximity measure, we conclude subsequential convergence to a stationary point.

Example 39 (Variable Euclidean Distance). A simple but far-reaching extension of Example 38 is the following. Let (A k ) k∈N be a sequence of symmetric positive definite matrices such that the smallest and largest eigenvalues are in [c 1 , c 2 ] for some 0 < c 1 < c 2 < +∞, i.e. 0 < inf k x, A k x < sup k x, A k x < +∞ for all x ∈ R N . Each matrix A k induces a metric on R N via the inner product x, Ax for x, x ∈ R N . The induced norm is a Bregman proximity function

D h k (x, x) = 1 2 |x -x| 2 A k := 1 2 x -x, A k (x -x) ,
generated analogously to Example 38. Except the boundedness of the eigenvalues of (A k ) k∈N there are no other restrictions. All the conditions mentioned in Example 38 are easily verified.

From now on, we restrict to iteration-independent Bregman distance functions, knowing that we can flexibly adapt the Bregman distance in each iteration.

Example 40 (Boltzmann-Shannon entropy). The Boltzmann-Shannon entropy is

D h (x, x) = N i=1 x (i) (log(x (i) ) -log(x (i) )) -(x (i) -x(i) )
where

x (i) denotes the i-th coordinate of x ∈ R N . D h is generated by the Legendre function h(x) = N i=1 x (i) log(x (i)
), which has the domain [0, +∞[ N . Since h is additively separable, w.l.o.g., we restrict the discussion to N = 1 in the following.

We verify Assumption 2. Let (x k ) k∈N and (x k ) k∈N be bounded sequences in int dom

h = ]0, +∞[ with x k -xk → 0 for k → ∞. For any convergent subsequence x k → x * as k K → ∞ for some K ⊂ N also xk → x * as k K → ∞ and x * ∈ [0, +∞[. Since h is continuous on cl dom h = [0, +∞[ (define h(0) = 0 log(0) = 0), D h (x k , xk
) → 0 for any convergent subsequence, hence for the full sequence. The same argument shows that the converse implication is also true, hence the Boltzmann-Shannon entropy satisfies Assumption 2.

For the model functions from Section 5.1, we show that Assumption 3 holds for x * ∈ int dom h, i.e. ∇h(x k ) -∇h(x k ) → 0 for sequence (x k ) k∈N and (x k ) k∈N with x k → x * and

x k -xk → 0 for k K → ∞ for some K ⊂ N. This condition is satisfied, because ∇h is continuous on int dom h, hence lim k K →∞ ∇h(x k ) = lim k K →∞ ∇h(x k ) = ∇h(x * ).
Since dom h = cl dom h, it suffices to verify Condition (iii) of Theorem 22 to guarantee subsequential convergence to a stationary point. For x * ∈ [0, +∞[ and a bounded sequence (ỹ k ) k∈N in int dom h as in Theorem 22, we need to show that D h (x * , ỹk ) → 0 as k K → ∞ for K ⊂ N such that ỹk → x * as k K → ∞. This result is clearly true for x * > 0, thanks to the continuity of log. For x * = 0, we observe x * log(ỹ k ) → 0 for k K → ∞, hence Condition (iii) of Theorem 22 holds, and subsequential convergence to a stationary point is guaranteed.

Example 41 (Burg's entropy). For optimization problems with non-negativity constraint, Burg's entropy is a powerful distance measure. The associated Bregman distance

D h (x, x) = N i=1 x (i) x(i) -log x (i) x(i) -1 is generated by the Legendre function h(x) = -N i=1 log(x (i)
) which is defined on the domain ]0, +∞[ N . Approaching 0, the function h grows towards +∞. In contrast to the Bregman functions in the examples above, Burg's entropy does not have a Lipschitz continuous gradient, and is therefore interesting for objective functions with the same deficiency.

W.l.o.g. we consider N = 1. Assumption 2 for two bounded sequences (x k ) k∈N and (x k ) k∈N in ]0, +∞[ reads

x k -xk → 0 ⇔ x k xk -log x k xk → 1 ,
which is satisfied if the limit points lie in ]0, +∞[ since

x k -xk → 0 ⇔ x k /x k → 1 for k K → ∞ and log is continuous at 1.
For the model functions in Section 5.1, Assumption 3 requires ∇h(x k ) -∇h(x k ) → 0 for sequence (x k ) k∈N and (x k ) k∈N in int dom h with x k → x * and x k -xk → 0 for k K → ∞ for some K ⊂ N. By continuity, this statement is true for any x * > 0. For x * = 0, the statement is in general not true. Also Condition (iv) in Theorem 22 can, in general, not be verified. Therefore, if a model functions is complemented with Burg's entropy, then the objective should be continuous on the cl dom h. Stationarity of limit points is obtained, if they lie in int dom h.

Applications

We discuss in this section some numerical experiments whose goal is to illustrate the wide applicability of our algorithmic framework. The applicability of our results follows from the considerations in Section 5. Actually, the considered objective functions are all continuous, i.e. Theorem 22(i) is satisfied.

Robust Non-linear Regression

We consider a simple non-smooth and non-convex robust regression problem [START_REF] Hampel | Robust Statistics: The Approach Based on Influence Functions[END_REF] of the form

F i (u) -y i 1 , F i (u) := P j=1 b j exp(-a j x i ) , (23)
where (x i , y i ) ∈ R × R, i = 1, . . . , M is a sequence of covariate-observation pairs. We assume that (x i , y i ) are related by y i = F i (u) + n i , where n i is the error term and u = (a, b) are the unknown parameters. We assume that the errors are iid with Laplacian distribution, in which case the data fidelity devised by a maximum likelihood argument is 1 -norm as used in [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF].

We define model functions by linearizing the inner functions F i as suggested by the model function in Example 35. Complemented by an Euclidean proximity measure (with τ > 0) the convex subproblem (5) to be solved inexactly is the following: ũ = argmin

u∈R P ×R P M i=1 K i u -y i 1 + 1 2τ |u -ū| 2 , y i := y i -F (ū) + K i ū ,
where K i := DF i (ū) : R P × R P → R is the Jacobian of F i at the current parameters ū. We solve the (convex) dual problem (cf. [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Combettes | Dualization of signal recovery problems[END_REF]) with warm starting up to absolute step difference 10 -3 . As mentioned in Remark 12, backtracking on τ could be used (cf. ProxDescent [START_REF] Lewis | A proximal method for composite minimization[END_REF]); denoted prox-linear and prox-linear2 in the following. This requires to solve the subproblem for each trial step. This is the bottleneck compared to evaluating the objective. The line search in Algorithm 2 only has to evaluate the objective value. This variant is denoted prox-linear-LS in the following. A representative convergence result in terms of the number of accumulated iterations of the subproblems is shown in Figure 1. For this random example, the maximal noise amplitude is 12.18, and the maximal absolute deviation of the solution from the ground truth is 0.53, which is reasonable for this noise level. Algorithm prox-linear-LS requires significantly fewer subproblem iterations than prox-linear and prox-linear2. For prox-linear2 the initial τ is chosen such that initially no backtracking is required.

For large scale problems, frequently solving the subproblems can be prohibitively expensive. Hence, ProxDescent cannot be applied, whereas our algorithm is still practical.

Image Deblurring under Poisson Noise

Let b ∈ R nx×ny represent a blurry image of size n x ×n y corrupted by Poisson noise. Recovering a clean image from b is an ill-posed inverse problem. It is a common problem, for example, in fluorescence microscopy and optical/infrared astronomy; see [START_REF] Bertero | Image deblurring with Poisson data: from cells to galaxies[END_REF] and references therein. A popular way to solve it is to formulate an optimization problem [START_REF] Zanella | Efficient gradient projection methods for edge-preserving removal of Poisson noise[END_REF] where A is a circular convolution (blur) operator. The first term (coined data term) in the objective f is the Kullback-Leibler divergence (Bregman distance generated by the Boltzmann-Shannon entropy x log(x)), which, neglecting additive constant terms, is given by f 1 (u) := D KL (b, Au) := i,j

(Au) i,j -b i,j log((Au) i,j ) , f 1 is well-suited for Poisson noise removal [START_REF] Vardi | A statistical model for positron emission tomography[END_REF]. The second term (coined regularization term) involves a penalty φ : R 2 → R applied to spatial finite differences (Du) i,j := ((Du) 1 i,j , (Du) 2 i,j ) in horizontal direction (Du) 1 i,j := u i+1,j -u i,j for all (i, j) with i < n x , and 0 otherwise; and vertical direction (Du) 2 i,j (defined analogously). The function φ in the regularization is usually chosen to favor "smooth" images with sharp edges. The relative importance of both the data and regularization terms is weighted by λ > 0.

For convex penalties φ, algorithms for solving problem [START_REF] Lee | Learning the part of objects from nonnegative matrix factorization[END_REF] are available (e.g. primaldual proximal splitting) provided that φ is simple (in the sense that its Euclidean proximal mapping can be computed easily). But if one would like to exploit the gradient of f 1 explicitly, things become more intricate. The difficulty comes from the lack of global Lipschitz continuity of ∇f 1 (u). A remedy is provided by Bauschke et al. [START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF]. They have shown that, instead of the global Lipschitz continuity, the key property is the convexity of Lh -f 1 for a Legendre function h and sufficiently large L, which can be achieved using Burg's entropy h(u) =i,j log(u i,j ) ([6, Lemma 7]).

However, non-convex penalties φ are known to yield a better solution [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF][START_REF] Blake | Visual Reconstruction[END_REF][START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF]. In this case, the algorithmic framework of Bauschke et al. [START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF] is not applicable anymore, whereas our framework is applicable. Due to the lack of strong convexity of Burg's entropy also the algorithm of Bonettini et al. [START_REF] Bonettini | Variable metric inexact line-search based methods for nonsmooth optimization[END_REF] cannot be used. Note that Burg's entropy is strongly convex on bounded subsets of ]0, +∞[, however, the subset cannot be determined a priori.

The abstract framework proposed in this paper appears to be the first algorithm with convergence guarantees for solving [START_REF] Lee | Learning the part of objects from nonnegative matrix factorization[END_REF] with a smooth non-convex regularizer.

In our framework, we choose φ : t ∈ R 2 → log(1 + ρ|t| 2 ), which is smooth but non-convex. The model functions are defined as in Example 33. We also use h as the Burg's entropy to generate the Bergman proximity function (see Example 41). Thus, the subproblems (5) which emerge from linearizing the objective f in (24) around the current iterate ū ũ = argmin

u∈R nx×ny u -ū, ∇f (ū) + 1 τ i,j
u i,j ūi,j -log u i,j ūi,j can be solved exactly in closed-form ũi,j = ūi,j /(1 + τ (∇f (ū)) i,j ūi,j ) for all i, j. A result for the successful Poisson noise removal and deblurring is shown in Figure 2.

Structured Matrix Factorization

Structured matrix factorization problems are crucial in data analysis. It has many applications in various areas including blind deconvolution in signal processing, clustering, source separation, dictionary learning, etc.. There is a large body of literature on the subject and we refer to e.g. [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-Way Data Analysis and Blind Source Separation[END_REF][START_REF] Chaudhuri | Blind Image Deconvolution[END_REF][START_REF] Starck | Sparse image and signal processing: wavelets, curvelets, morphological diversity[END_REF][START_REF] Xu | A survey of dictionary learning algorithms for face recognition[END_REF] and references therein for a comprehensive account.

The problem. Given a data matrix A ∈ R M ×N whose N M -dimensional columns are the data vectors. The goal is to find two matrices U ∈ R M ×K and Z ∈ R K×N such that

A = U Z + Q ,
where Q ∈ R M ×N accounts for an unknown error. The matrices U and Z (called also factors) enjoy features arising in a specific application at hand (see more below).

To solve the matrix factorization problem, we adopt the optimization approach and we consider the non-convex and non-smooth minimization problem min

U ∈U ,Z∈Z f (A, U Z) + λg(Z) , f (A, U Z) := 1 2 A -U Z 2 F . (25)
The term f (A, U Z) stands for proximity function that measures fidelity of the approximation of A by the product U Z of the two factors. We here focus on the classical case where the fidelity is measured via the Frobenius norm • F , but other data fidelity measures can also be used just as well in our framework, such as divergences (see [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-Way Data Analysis and Blind Source Separation[END_REF] and references therein). The sets U, Z, which are non-empty closed and convex, and the function g ∈ Γ 0 are used to capture specific features of the matrices U and Z arising in a specific application as we will exemplify shortly. The influence of g is weighted by the parameter λ > 0. Many (if not most) algorithms to solve the matrix factorization problem (25) are based on Gauss-Seidel alternating minimization with limited convergence guarantees [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-Way Data Analysis and Blind Source Separation[END_REF][START_REF] Chaudhuri | Blind Image Deconvolution[END_REF][START_REF] Starck | Sparse image and signal processing: wavelets, curvelets, morphological diversity[END_REF] 4 . The PALM algorithm proposed recently by Bolte et al. [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF], was designed specifically for the structure of the optimization problem [START_REF] Lewis | A proximal method for composite minimization[END_REF]. It can then be successfully applied to solve instances of such a problem with provably guaranteed convergence under some assumptions including the Kurdyka-Lojasiewicz property. However, though it can handle non-convex constraint sets and functions g, it does not allow to incorporate Bregman proximity functions.

In the following, we show how our algorithmic framework can be applied to a broad class of matrix factorization instances. In particular, a distinctive feature of our algorithm is that it can naturally and readily accommodate for different Bregman proximity functions and it has no restrictions on the choice of the step size parameters (except positivity). A descent is enforced in the line search step, which follows the proximal step.

A generic algorithm. We apply Algorithm 1 to solve this problem, where the model functions are chosen to linearize the data fidelity function f (A, U Z), according to Example 33. The convex subproblems to be solved in the algorithm have the following form:

( Ũ , Z) = argmin U ∈U ,Z∈Z λg(Z)+ Z -Z, Ū ( Ū Z -A) F + D h Z (Z, Z) + U -Ū , ( Ū Z -A) Z F + D h U (U, Ū )
where •, • F stands for the Frobenius inner product. The Bregman proximity functions D h Z (•, •) and D h U (•, •) provide the flexibility to handle a variety of constraint sets U and Z. In the following, we list different choices for the constraint sets and explain how to incorporate them into the optimization procedure. Due to the structure of the optimization problem, the variables U and Z can be handled separately. The only coupling is the data fidelity function f , which is linearized and therefore easy to incorporate.

Examples of constraints U. There are many possible choices for the set U depending on the application at hand.

• Unconstrained case:

U 1 = R M ×K .
In the unconstrained case, a suitable Bregman proximity function is given by the Euclidean distance D h U (U, Ū ) = 1 2τ U U -Ū 2 F with step size parameter τ U . The resulting update step with respect to the dictionary U is a gradient descent step.

• Zero-mean and normalization:

U 2 = U ∈ R M ×K | ∀j : M i=1 U 2 i,j ≤ 1 , ∀j ≥ 2 : M i=1 U i,j = 0 .
This choice of the constraint set leads to a natural normalization of the columns of U that removes the scale ambiguity due to bilinearity. This choice is very classical in dictionary learning, see, e.g., [START_REF] Xu | A survey of dictionary learning algorithms for face recognition[END_REF]. As in dictionary learning, the average of the first column may not be enforced to be zero, in order to allow the first column to absorb the mean value of the data points.

By separability of U 2 , the Euclidean projection onto it is simple. This projector is column-wise achieved by subtracting the mean, and then projecting the result onto the Euclidean unit ball. Thus we advocate D h U (U, Ū ) = 1 2τ U U -Ū 2 F with step size parameter τ U . In turn, the subproblem with respect to U amounts to a projected gradient descent step.

• Non-negativity and normalization:

U 3 = U ∈ R M ×K | ∀j : M i=1 U i,j = 1 , ∀i, j : U i,j ≥ 0 .
This choice is adopted in non-negative matrix factorization (NMF) [START_REF] Lee | Learning the part of objects from nonnegative matrix factorization[END_REF]. The constraint set U 3 is column-wise a unit simplex constraint. This constraint can be conveniently handled by choosing D h U (U, Ū ) = 1 τ U i,j U i,j (log(U i,j ) -log( Ūi,j )) -U i,j + Ūi,j , which is the Bregman function generated by the entropy h U (U ) = 1 τ U i,j U i,j log(U i,j ) with step size parameter τ U . This is a more natural choice than the Euclidean proximity distance. Indeed, the update step with respect to U results in Ũi,j = Ūi,j exp(-τ U (C U ) i,j ) M p=1 Ūp,j exp(-τ U (C U ) p,j ) ∀i = 1, . . . , M ; ∀j = 1, . . . , K ,

where we use the shorthand notation C U := ∇ U f (A, Ū Z) = Ū ( Ū Z -A) for the partial gradient of f with respect to U . The exponential function is applied entry-wise, hence naturally preserving positivity. Note that the Euclidean projector onto U 3 necessitates to compute the projector on the simplex which can be achieved with sorting [START_REF] Michelot | A finite algorithm for finding the projection of a point onto the canonical simplex of R n[END_REF].

Examples of constraints Z. There are also several possible choices for the set Z and regularizing function g depending on the application at hand.

• Unconstrained case: Z 1 = R K×N and g(Z) = 0 .

This case can be handled using a gradient descent step, analogously to the related update step with the constraint set U 1 .

• Non-negativity:

Z 2 = Z ∈ R K×N | ∀i, j : Z i,j ≥ 0 and g(Z) = 0 .
This constraint is used in conjunction with U 3 in NMF. It can be handled either with a Euclidean proximity function (which amounts to projecting on the non-negative orthant), or via a Bregman proximity function D h Z (Z, Z) generated by the Boltzmann-Shannon entropy (h Z (Z) = 1 τ Z i,j Z i,j log(Z i,j )) or, alternatively, Burg's entropy (h Z (Z) = -1 τ Z i,j log(Z i,j )), with step size parameter τ Z . The update with respect to Z then reads Zi,j = Zi,j exp(-τ Z (C Z ) i,j ) ∀i = 1, . . . , K ; ∀j = 1, . . . , N , where we use the shorthand notation C Z := ∇ Z f (A, Ū Z) = ( Ū Z -A) Z for the partial gradient of f with respect to Z.

• Sparsity constraints:

Z 3 = R K×N and g(Z) = Z 1 .
The introduction of sparsity has been of prominent importance in several matrix factorization problems, including dictionary learning [START_REF] Olshausen | Sparse coding with an overcomplete basis set: A strategy employed by V1?[END_REF], NMF [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF] 5 and source separation [START_REF] Starck | Sparse image and signal processing: wavelets, curvelets, morphological diversity[END_REF]. The Euclidean proximal mapping of the 1 -norm is the entry-wise softthresholding, hence giving the update step with respect to Z as Zi,j = max{0, 1-λτ Z /| Zi,j -τ Z (C Z ) i,j |}( Zi,j -τ Z (C Z ) i,j ) ∀i = 1, . . . , K ; ∀j = 1, . . . , N .

• Low rank constraint:

Z 3 = R K×N and g(Z) = Z * .
The nuclear norm or 1-Schatten norm Z * is the sum of the singular values. It is known to be the tightest convex relaxation to the rank and was shown to promote low rank solutions [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF]. Such a regularization would be useful in the situation where columns of A are (to a good approximation) clustered on a few linear subspaces spanned by the columns of U , i.e. the columns of A can be explained by columns of U from the same subspace ("cluster").

The Euclidian proximal mapping of the nuclear norm is the soft-thresholding applied to the singular values. In turn, the update step with respect to Z reads Zi,j = W diag((max{0, 1 -λτ Z /σ i }σ i ) i )V , where W , V are respectively the matrices of left and right singular vectors of Z -τ Z C Z , and σ is the associated vector of singular values.

We have presented an algorithmic framework, that unifies the analysis of several first order optimization algorithms in non-smooth non-convex optimization such as Gradient Descent, Forward-Backward Splitting, ProxDescent, and many more. The algorithm combines sequential Bregman proximal minimization of model functions, which is the key concept for the unification, with an Armijo-like line search strategy. The framework reduces the difference between algorithms to the model approximation error measured by a growth function. For the developed abstract algorithmic framework, we establish subsequential convergence to a stationary point and demonstrate its flexible applicability in several difficult inverse problems from machine learning, signal and image processing.

Remark 21 .

 21 Assumption 3 is common for abstract algorithms. Attouch et al. [2], for example, use a relative error condition of the form |∇f |(x k+1 ) ≤ b|x k+1 -x k |, b ∈ R. A weaker sufficient condition for Assumption 3 is |∇f |(x k+1 ) ≤ ψ(|x k+1 -x k |) for some continuous function ψ : R + → R + with ψ(0) = 0; See Corollary 25. See also Remark 24. For explicit examples, we refer to Section 5.1.

  and by Assumption 3, the slope vanishes, hence the limiting slope |∇f |(x * ) at x * also vanishes, which concludes the proof.

  min u:=(a,b)∈R P ×R P M i=1

Figure 1 :

 1 Figure 1: Objective value vs. accumulated number of subproblem iterations for (23).

Figure 2 :

 2 Figure 2: Deblurring and Poisson noise removal by solving (24). From left to right: clean, noisy, and reconstructed image (PSNR: 25.86).

It is often easy to find a feasible point. Of course, there are cases, where finding an initialization is a problem itself. We assume that the user provides a feasible initial point.

Note that inf k η k > 0 is equivalent to lim inf k η k > 0 as we assume η k > 0 for all k ∈ N.

The example is not meant to be meaningful and the model function to be algorithmically the best choice. This example shall demonstrate the flexibility and problem adaptivity of our framework.

For very specific instances, a recent line of research proposes to lift the problem to the space of low-rank matrices, and then use convex relaxation and computationally intensive conic programming that are only applicable to small-dimensional problems; see, e.g.,[START_REF] Ahmed | Blind deconvolution using convex programming[END_REF] for blind deconvolution.

Strictly speaking, Z 3 should be the non-negative orthant for sparse NMF. But this does not change anything to our discussion since computing the Euclidean proximal mapping of the 1 norm restricted to the non-negative orthant is easy.