Bregman Itoh–Abe Methods for Sparse Optimisation

In this paper we propose optimisation methods for variational regularisation problems based on discretising the inverse scale space flow with discrete gradient methods. Inverse scale space flow generalises gradient flows by incorporating a generalised Bregman distance as the underlying metric. Its discrete-time counterparts, Bregman iterations and linearised Bregman iterations are popular regularisation schemes for inverse problems that incorporate a priori information without loss of contrast. Discrete gradient methods are tools from geometric numerical integration for preserving energy dissipation of dissipative differential systems. The resultant Bregman discrete gradient methods are unconditionally dissipative and achieve rapid convergence rates by exploiting structures of the problem such as sparsity. Building on previous work on discrete gradients for non-smooth, non-convex optimisation, we prove convergence guarantees for these methods in a Clarke subdifferential framework. Numerical results for convex and non-convex examples are presented.


Introduction
We consider the constrained optimisation problem min x∈C V (x), (1.1) for an objective function V : R n → R and constraint C ⊂ R n .The function V may be non-convex and non-smooth, as outlined in Assumption 3.1.In this paper, we propose and study optimisation schemes by using tools from geometric numerical integration to solve the inverse scale space (ISS) flow.
The ISS flow is a differential system which generalises gradient flows by replacing the Euclidean distance by a Bregman distance, defined via a convex Bregman (distance generating) function J : R n → R. The ISS flow is given by ṗ(t) = −∇V (x(t)), p(t) ∈ ∂ J (x(t)). (1.2) Here, ∂ J is the convex subdifferential of J , to be defined in Sect.2, and R := R ∪ {±∞}.The term inverse scale space flow goes back to Scherzer and Groetsch [50].This flow is typically derived as the continuous-time limit of Bregman iterations-methods for solving variational regularisation problems.Like gradient flows, the ISS flow is a dissipative system, and its dissipative structure is determined by the function J .This allows one to solve (1.1) while incorporating a priori information into the optimisation scheme, with the benefits of converging to superior solutions, and doing so faster.For background on the ISS flow and Bregman distances, see Sect.2.1.Geometric numerical integration deals with numerical integration methods that preserve geometric structures of the continuous system.Such geometric structures include dissipation or conservation of energy and Lyapunov functions.In recent years, geometric numerical integration-and numerical integration in general-has gained interest within the mathematical optimisation community, as a framework for formulating iterative schemes that are dissipative or amenable to Lyapunov arguments.
We propose to discretise the inverse scale space flow with discrete gradient methods.These are methods from geometric numerical integration that preserve the aforementioned geometric structures in a general setting.In recent papers [21,26,45,46], optimisation schemes based on discretising gradient flows with discrete gradients have been analysed and implemented for various problems.Favourable properties of the discrete gradient methods include unconditional dissipation, i.e. dissipation is ensured for any time step, and the Itoh-Abe discrete gradient method is derivative-free and has convergence guarantees in the non-smooth, non-convex setting [45].For smooth problems, the theoretical convergence rates of the discrete gradient methods match those of explicit gradient descent and coordinate descent [21].The drawback of these methods is that the updates are in general implicit.Nevertheless, for many simple variational problems, the updates turn out to be explicit.
In this paper, we study the Itoh-Abe discrete gradient method applied to the ISS flow.We prove that the method is well defined and converges to a set of stationary points for non-smooth, non-convex functions.Furthermore, building on the paper by Miyatake et al. [38] which pointed out the equivalence between the discrete gradient methods for linear systems and successive-over-relaxation (SOR) methods, we point out equivalencies of various approaches to least squares problems.
Bregman iterations, and related methods, are closely tied to inverse problems and regularisation methods, particularly in signal processing.We consider numerical examples in this setting as well.

Related Literature
Spurred by applications for variational regularisation in image processing and compressed sensing, the ISS flow and the Bregman method have been active areas of research during the last decade.The Bregman iterative method was originally proposed by Osher et al. [41] in 2005 for total variation-based image denoising, representing an extension of the Bregman proximal algorithm [14,19,32,55] to nonsmooth Bregman functions.Subsequently the ISS flow was derived and analysed by Burger et al. [7,[9][10][11].Since then, researchers have studied the ISS flow with applications to generalised spectral analysis in a nonlinear setting, i.e. by Burger et al. [8], Gilboa et al. [23], and Schmidt et al. [51].The Bregman method has been studied for 1 -regularisation and compressed sensing by Goldstein and Osher [24] and Yin et al. [59], and extended to primal-dual algorithms by Zhang et al. [61].
The linearised Bregman method was proposed by Yin et al. [59] for applications to 1 -regularisation and compressed sensing, and further studied in this setting by Cai et al. [12], and Dong et al. [18].An extension for non-convex problems was proposed by Benning et al. [3], proving global convergence for functions that satisfy the Kurdyka-Łojasiewicz property.Lorenz et al. [34,52] proposed a sparse variant of the Kaczmarz method for linear problems based on linearised Bregman iterations.These and other methods were unified in a Split Feasibility Problems framework for general convergence results by Lorenz et al. [33].For further details on Bregman iterations and linearised Bregman methods, we refer to [4].
We review papers on discrete gradient methods for optimisation based on gradient flows.Grimm et al. [26] used discrete gradients to solve variational regularisation problems in image analysis and proved convergence to a set of stationary points for continuously differentiable functions.Ehrhardt et al. [21] provided additional analysis for the methods, including convergence rates for smooth, convex problems and Polyak-Łojasiewicz functions, as well as well-posedness of the implicit equation.Ringholm et al. [46] applied the Itoh-Abe discrete gradient method to non-convex image problems with Euler's elastica regularisation.
Furthermore, several recent works have looked at discrete gradient methods in more general settings.Riis et al. [45] studied the Itoh-Abe discrete gradient method in the setting of derivative-free optimisation of non-smooth, non-convex objective functions, and proved that the method converges to a set of stationary points in the Clarke subdifferential framework.Celledoni et al. [13] extended the Itoh-Abe discrete gradient method to optimisation problems defined on Riemannian manifolds.Hernández-Solano et al. [29] combined a discrete gradient method with Hopfield networks in order to preserve a Lyapunov function for optimisation problems.
We point out that a central feature of discrete gradient methods is that due to their implicit formulation, no restrictions are required for the time steps.This will also be the case for the analysis in this paper.
Beyond discrete gradients, other methods from numerical integration have led to notable developments in optimisation in recent years.Recent papers by Su et al. [54] and Wibisono et al. [56] study second-order ODEs based on the continuoustime limit of Nesterov's accelerated gradient descent [39], in order to gain insight into the acceleration phenomenon.In this setting, Wilson et al. [57] provided a framework for Lyapunov analysis of optimisation schemes and continuous-time dynamics, and Betancourt et al. [5] proposed a framework for symplectic integration for optimisation.On a related note, Scieur et al. [53] showed that several accelerated optimisation schemes can be derived as multi-step integration schemes from numerical analysis.An optimisation scheme based on numerically integrating dissipative Hamiltonian conformal systems was proposed by Maddison et al. [35], with the aim of ensuring linear convergence for a larger group of functions than classical Other numerical integration methods of relevance include implicit Runge-Kutta methods, where energy dissipation is ensured under mild time step restrictions [27], and explicit stabilised methods for solving strongly convex problems [20].Another example is the study of gradient flows in metric spaces and minimising movement schemes [1], concerning gradient flow trajectories under other measures of distance, such as the Wasserstein metric [49].

Structure and Contributions
The paper is structured as follows.For the remainder of Sect. 1, we define mathematical notation.In Sect.2, we provide preliminary material for convex and non-convex analysis, and introduce Bregman distances.In Sect.3, we propose a Bregman discrete gradient method based on the ISS flow, and prove well-posedness and convergence results in a non-convex, non-smooth framework.In Sect.4, we discuss particular examples of Bregman discrete gradient methods.In Sect.5, we present results from numerical experiments.In Sect.6, we conclude.

Notation and Preliminaries
Throughout the paper, we denote by • the 2 -norm, i.e.
• p denotes the usual p -norm.For ε > 0 and x ∈ R n , we denote by B ε (x) the open ball {y ∈ R n : y − x < ε}.We denote by {e 1 , . . ., e n } the standard coordinate vectors in R n .We define the extended reals as R := R ∪ {±∞}.
For a sequence (x k ) k∈N ⊂ R n , we denote by S its limit set, i.e. the set of accumulation points, For a matrix A ∈ R m×n , we denote by a i ∈ R n its ith row, and A T its transpose.Accordingly, a i j refers to the jth element of the ith column of A.
For x ∈ R, the sign operator sgn : R → {±1, 0} is defined as

Preliminaries for Convex and Non-convex Analysis
In this section, we provide preliminary material on nonsmooth analysis.In Sect.2.1, we cover convex analysis and Bregman distances, while in Sect.2.2 we cover nonconvex, nonsmooth analysis in the Clarke subdifferential framework.

Convex Analysis
We consider functions J : R n → R that are convex, proper, and lower-semicontinuous (see [22,48] for details on this class of functions).

Definition 2.1 (Effective domain)
The effective domain of a function J : R n → R is defined as dom(J ) = {x ∈ R n : J (x) < ∞}.

Definition 2.2 (Subgradients and subdifferentials)
The subdifferential of a convex function J : R n → R at x ∈ R n is the set of vectors Vectors in ∂ J (x) are called subgradients of J at x.
For the theory of convex functions and their subdifferentials, see [48].
We consider the characteristic function χ C of a convex, closed set C ⊂ R n , defined by This is a convex, proper, lower-semicontinuous function, and We can now define the generalised Bregman distance [6] of a convex function J .While Bregman distances are non-negative due to convexity of J , they are not metrics as they do not generally satisfy symmetry or a triangle inequality.However, a related object is the symmetric Bregman distance.Definition 2.5 (Symmetric Bregman distance) Given p ∈ ∂ J (x) and q ∈ ∂ J (y), the symmetric Bregman distance between x and y is given by If μ > 0, J is said to be strongly convex.
For strongly convex Bregman functions, the following property of Bregman distances is immediate.

Non-Convex Subdifferential Analysis
We summarise the main concepts of the Clarke subdifferential for locally Lipschitz continuous, non-smooth, nonconvex functions V : R n → R, and refer to [16] for further details.
Definition 2.8 (Lipschitz continuity) V is Lipschitz of rank L near x if there exists ε > 0 such that for all y, z ∈ B ε (x), one has V is locally Lipschitz continuous if the above property holds for all x ∈ R n .
Definition 2.9 (Clarke directional derivative) For V Lipschitz near x and for a vector d ∈ R n , the Clarke directional derivative is given by Definition 2.10 (Clarke subdifferential) Let V be locally Lipschitz continuous and x ∈ R n .The Clarke subdifferential of V at x is given by The Clarke subdifferential was introduced by Clarke in [15].
In this paper, we will consider constrained optimisation problems arg min x∈C V (x) for a closed, convex subset C. We assume in the convergence analysis that C is box constraints of the form To have a formal notion of first-order optimality in this setting, we define the tangent cone.
As the following proposition shows, the above definition generalises minimisers for convex functions, and Clarke stationary points in the unrestricted case.

The Discrete Gradient Method for the Inverse Scale Space Flow
In what follows, we discuss the inverse scale space flow, and propose a Bregman discrete gradient method by using discrete gradients to solve the differential system.

Inverse Scale Space Flow and Bregman Methods
For a convex function J : R n → R, objective function V : R n → R and starting points x(0) = x 0 ∈ R n , p(0) ∈ ∂ J (x 0 ), the ISS flow is the dissipative differential system given by (1.2).If J were twice continuously differentiable and μ-convex, then (1.2) could be rewritten as and the energy V (x(t)) would dissipate over time as We briefly discuss variants of Bregman methods as discretisations of (1.2).The Bregman method is derived by backward Euler discretisation of (1.2) and is given by which can be rewritten as From (3.1), we see that the Bregman method is dissipative, as Similarly, the linearised Bregman method is derived by forward Euler discretisation of (1.2) and is given by or equivalently The ISS flow and Bregman methods are considered for solving ill-conditioned linear systems Ax = b, with the objective function In this setting, iterates of both the Bregman method and the linearised Bregman method converge [4,33] to a solution of Furthermore, applications of the ISS flow include image denoising with reduced contrast-loss and staircasing effects [41], recovering eigenfunctions [51], and identifying sparse or low-rank structures [59].
We make the following assumptions on the objective function V , the constraints C, as well as the Bregman function J .C if and only if for all coordinate vectors e i , we have and μ-convex with μ > 0. Furthermore, , where j i : R → R are convex and Lipschitz continuous on [l i , u i ] for each i.

The Bregman Discrete Gradient Method
In what follows, we define discrete gradients and propose the Bregman discrete gradient method by discretising the ISS flow.
Discrete gradients are tools from geometric numerical integration.In geometric integration, one studies methods for numerically solving ODEs while preserving certain structures of the continuous system-see [28,36] for an introduction.Discrete gradients are tools for solving first-order ODEs that preserve energy conservation laws, dissipation laws, and Lyapunov functions [25,30,37,44].
The Itoh-Abe discrete gradient [30] (also known as the coordinate increment discrete gradient) is given by where 0/0 is interpreted as The Itoh-Abe discrete gradient is derivative-free and is evaluated by successively computing difference quotients.For optimisation problems where the gradient is expensive to compute (e.g.[46]), or unavailable (e.g.[45]), Itoh-Abe discrete gradient methods are useful for ensuring convergence guarantees and competitive convergence rates [21] without requiring derivatives.
We propose the Bregman discrete gradient method as follows.For a starting point x 0 ∈ R n , subgradient p 0 ∈ ∂ J (x 0 ), and time steps (τ k ) i∈N > 0, solve for k = 0, 1, . .., This scheme preserves the dissipative structure of the ISS flow and (linearised) Bregman methods, as we see by applying the mean value property (3.2) and (3.5). V Furthermore, if we plug in J (x) = x 2 /2, we recover the discrete gradient method for gradient flows [21].Observe that the dissipative structure holds for arbitrary positive time steps τ k > 0. We thus only require that τ k ∈ [τ min , τ max ] for arbitrary bounds τ min , τ max > 0.
By Assumption 3.1 d), the subdifferential of J is separable in the coordinates, i.e.
It follows that solving the Bregman Itoh-Abe equation (3.5) corresponds to successively solving n scalar inclusions, i.e. given x k , and p k , we want x k+1 and p k+1 that solve Remark 3.3 While discrete gradients are not defined for nonsmooth functions V , the resultant Bregman discrete gradient method (3.7) is still well defined, provided the function V has a well defined Clarke subdifferential.This is the case, e.g. if V is locally Lipschitz continuous [16].Furthermore, note that the properties used to derive the dissipative structure (3.6) holds in this setting.

Theoretical Results
In this section, we prove that a solution to the Bregman discrete gradient method (3.7) exists, and provide conditions in which the update is unique.Furthermore, we prove that the accumulation points of the iterates (x k ) k∈N from the Bregman Itoh-Abe method are Clarke stationary.

Well-Posedness and Preliminary Results
The Bregman Itoh-Abe scheme (3.7) is in general implicit.We therefore want to ensure that an update exists and is easy to compute.Furthermore, we want to know under what conditions the updates are unique.In what follows, we address these questions.The corresponding analysis generalises that from the previous papers on discrete gradient methods for optimisation [21,45].
We first show that an update (not necessarily unique) always exists.Lemma 4.1 (Existence) For any τ > 0, x k ∈ R n , and p k ∈ ∂ J (x k ), there exists an update (x k+1 , p k+1 ) that satisfies (3.7).
Proof As (3.7) consists of successive scalar updates, it is sufficient to consider a scalar problem, v : R → R, j : R → R. For x ∈ R and p ∈ ∂ j(x) we either want y = x such that or y = x and p − τ w ∈ ∂ j(x), for some w ∈ ∂v(x).
If such a w exists, we are done.Otherwise, we have min{v o (x; 1), v o (x; −1)} < 0 and may assume that v o (x; 1) < 0. In this case, we will show that there exists y > x such that (4.1) holds.
Since p − τ v o (x; 1) > p and p ∈ ∂ j(x), we deduce that p − τ v o (x; 1) > ∂ j(x).By the outer semicontinuity of subdifferentials and definition of Clarke directional derivatives, there is δ > 0 such that On the other hand, as v is bounded below, is bounded below for all y ∈ [x + δ, +∞), while by μconvexity of j, we have ∂ j(y) ≥ ∂ j(x) + μ(y − x) for all y ∈ [x + δ, +∞).Hence, there is r 0 such that By continuity of v, and by outer semicontinuity of subdifferentials, it follows that there exists y ∈ (x + δ, x + r ) that solves (4.1).This concludes the proof.
Next we give conditions in which the update is unique.While this holds for any time step if the objective function is convex, it may fail for nonconvex functions, as discussed in [21, Section 5].
Proof Case (i) Suppose V is convex.The existence of a solution to (3.7) is guaranteed by Lemma 4.1.To establish uniqueness, we argue as follows.An update y k,i must satisfy The left-hand side is non-increasing with respect to x k+1 i due to the difference quotient term of the convex function V , while the right-hand side is strictly increasing due to the strong convexity of J .Hence there cannot be two distinct solutions for y k,i i to the scalar equation.This implies uniqueness of the update.
Case (ii): Suppose V (•) + η • 2 is convex and τ < μ/η.We rewrite (4.2) as Since τ η < μ, the function J (•) − τ η • 2 /2 is strongly convex, and therefore the right-hand side is strictly increasing with respect to x k+1 i .On the other side, as i is non-decreasing with respect to x k+1 i .Therefore, by the same reasoning as in the first case, the update to the Bregman Itoh-Abe method must be unique.
Note that while the convexity criteria in Proposition 4.2 are global, they are often generalised to local properties such as lower C 2 -smoothness [47] and prox-regularity [43].However, an analysis of these properties in the context of Itoh-Abe discrete gradient methods is beyond the theoretical and practical scope of this paper.We also add that lack of uniqueness is not an issue for the implementation or computational tractability of the scheme.
The following lemma summarises several useful properties of the iterates (x k , p k ) k∈N , and generalises Lemmas 3.3, 3.4 in [45].
then there exists a convergent subsequence of (x k , p k ) k∈N .(iv) The set of limit points S is compact, connected, and has empty interior.Furthermore, V is single-valued on S.
Proof Property (i) follows from (3.6).As V is bounded below and (V (x k )) k∈N is decreasing, V (x k ) → V * for some limit V * .Therefore, by (3.6),

Amended Scheme
In the next subsection, we prove that accumulation points of the iterates from the Bregman Itoh-Abe method are Clarke stationary.However, we first modify the Bregman Itoh-Abe method to 'forget' subgradients induced by the constraints.
We define the modified scheme as follows.For a starting point x 0 , p 0 ∈ ∂ J (x 0 ), and k ∈ N, update Observe that since 0 ∈ ∂χ C (x) for all x ∈ C, we still have p k+1 ∈ ∂ J (x k+1 ).It is also straightforward to verify that the previous results and analysis of dissipative structure in this section also hold for (4.3).
The reason for introducing the modified scheme (4.3) is so that the subgradient iterates p k i do not grow arbitrarily large when x k i equals u i or l i for several updates (recall that when the constraint is active, ∂χ [l i ,u i ] is unbounded).Otherwise, the iterates might get stuck at a non-stationary point for several iterations, since x k i is unable to leave {l i , u i } until the subgradient update in ∂χ [l i ,u i ] vanishes, yielding inefficient progress.Furthermore, this leads to pathological, albeit unlikely, examples where the iterates of (3.7) converge to non-stationary accumulation points.For completeness, we give such an example in Section A.

Main Convergence Theorem
Having introduced a modified Bregman Itoh-Abe scheme in the previous section, we proceed to state and prove the main theorem of this paper, namely that all accumulation points of the scheme (4.3) are non-stationary.We note that this also holds for the original Bregman Itoh-Abe method if the iterates (x k ) k∈N converge to a unique limit.Proof Let x * ∈ S and consider a convergent subsequence (x k j ) j∈N .We want to show for each basis vector e i that either V o (x * ; e i ) ≥ 0 or x * i = u i , and analogously that either V o (x * ; −e i ) ≥ 0 or x * i = l i .As the arguments are equivalent, we only prove the first case.
Suppose for contradiction that V o (x * ; e i ) < −η for some η > 0, and that x * i < u i .By the definition of the Clarke directional derivative, Definition 2.9, there are ε, δ > 0 such that for all x ∈ B ε (x * ) and λ ∈ (0, δ), we have Since x k j → x * and x k j +1 − x k j → 0, for each N ∈ N there exists K such that for all j ≥ K , we have and arrive at a contradiction.Thus, x * is a Clarke stationary point restricted to C.

Examples of Bregman Itoh-Abe Discrete Gradient Schemes
In this section, we describe several schemes based on the Bregman Itoh-Abe discrete gradient scheme (3.7).We will primarily consider objective functions of the form where A ∈ R n×n is a symmetric, positive semi-definite matrix, with strictly positive entries a i i > 0 on the diagonal.We are particularly interested in problems with underlying sparsity and/or constraints, with applications in image analysis.Throughout this section, we use a time step vector τ k coordinate-wise scaled by the diagonal of A, i.e.
for all k ∈ N, and some τ > 0.
We first introduce some well-known coordinate descent schemes for solving linear systems, which Miyatake et al. [38] showed were equivalent to the Itoh-Abe discrete gradient method.The successive-over-relaxation (SOR) method [60] updates each coordinate sequentially according to the rule where ω ∈ (0, 2).For ω = 1, this is the Gauss-Seidel method [60].The SOR method is equivalent to the Itoh-Abe discrete gradient method with V given by (5.1) with the time steps

Sparse SOR Method
We consider underdetermined linear systems and want to find sparse solutions x * .Hence we seek to apply the Bregman Itoh-Abe method (3.7) with objective function V given by (5.1), and for γ > 0. We term this the Bregman SOR (BSOR) method.By Proposition 4.2, the updates of this method are well defined and unique.One can verify that the updates are given as follows.Denote by xk+1 i the standard SOR coordinate update from x k i , (5.2).Furthermore, for p k ∈ ∂ J (x k ), we write ( Here S : R × R → R denotes the shrinkage operator S(x, λ) = sgn(x) max{|x| − λ, 0}, applied elementwise to a vector x.

Sparse, Regularised SOR
If b = Ax true + , where x true is the sparse ground truth and is noise, then it may be necessary to regularise the objective function as well.Hence we consider the objective function for some regularisation parameter λ > 0. The nonsmoothness induced by • 1 satisfies Assumption 3.1, so Theorem 4.5 implies Bregman Itoh-Abe discrete gradient methods will still converge to stationary points of this problem.For both J (x) = 1 2 x 2 + γ x 1 and J (x) = 1 2 x 2 , the scheme (3.7) can be expressed in closed form for (5.5), on a case-by-case basis.However, for purposes of brevity, we leave including this in "Appendix".

Equivalence of Iterative Methods for Linear Systems
In what follows, we discuss and demonstrate equivalencies for different iterative methods for solving linear systems.We recall from the previous section that the SOR method (5.2) is equivalent to the Itoh-Abe discrete gradient method [38].
The explicit coordinate descent method [2,58] for minimising V is given by where α i > 0 is the time step.As mentioned in [58], the SOR method is also equivalent to the coordinate descent method with V given by (5.1) and the time step α i = ω/a i i .Hence, in this setting, the Itoh-Abe discrete gradient method is equivalent not only to SOR methods, but to explicit coordinate descent.
It is not surprising that these iterative coordinate methods turn out to be the same, given that the gradient V in (5.1) is linear.Furthermore, these equivalencies extend to discretisations of the inverse scale space flow with J given by (5.3).The resultant Bregman Itoh-Abe scheme for (5.1) is described in Sect.5.1.We may compare this to a Bregman linearised coordinate descent scheme, One can verify that this scheme is equivalent to (6.1) for the parameters

Numerical Examples
In this section, we present numerical results for the schemes described in Sect. 5.

Sparse SOR
We construct a matrix A ∈ R 1024×1024 from independent standard (zero mean, unit variance) Gaussian draws, and construct the sparse ground truth x true by choosing 10% of the indices at random determined by uniform draws on the unit interval.We then solve the problem arg min where b = Ax true .We compare the SOR method (J (x) = x 2 /2) and the BSOR method (J (x) = x 2 /2 + γ x 1 ), where γ = 1.We set time steps to τ = 2/diag (A), corresponding to the Gauss-Seidel method.See Fig. 1 for the results.
For the same test problem, but where the ground truth is binary, i.e.only takes values 1 or 0, see Fig. 2.

Sparse, Regularised SOR
We construct A ∈ R 1024×1024 and x true as in the previous subsection.However, we add noise to the data, i.e. b = Ax true + , where is independent Gaussian noise with a Fig. 3 Comparison of SOR and sparse SOR methods, for Gaussian linear system with noise.Top: Convergence rate for relative objective.Bottom: Support error with respect to iterates standard deviation of 0.1 Ax true ∞ .Since the added noise destroys the sparsity structure of A −1 b, the sparse SOR method fails to improve the convergence rate.The results for V (x) = Ax − b 2 /2 are in Fig. 3.
We therefore include regularisation in the objective function of the form where λ = 100, and with initialisation x 0 constructed by random, independent Gaussian draws.The results are visualised in Fig. 4.

Student-t Regularised Image Denoising
We consider a nonconvex image denoising model, previously presented in [40], given by Here {K i } N i=1 is a collection of linear filters, (ϕ i ) N i=1 ⊂ [0, ∞) are coefficients, Φ : R n → R is the nonconvex function based on the student-t distribution, defined as As impulse noise only affects a small subset of pixels, we use the data fidelity term x → x − x δ 1 to promote sparsity of x * −x δ for x * ∈ arg min F(x).As linear filters, we consider the simple case of finite difference approximations to firstorder derivatives of x.We note that by applying a gradient flow to this regularisation function, we observe a similarity to Perona-Malik diffusion [42].
We consider a Bregman Itoh-Abe method (abbreviated to BIA) for solving min x∈R n F(x) with the Bregman function to account for the sparsity of the residual x * −x δ , and compare the method to the regular Itoh-Abe discrete gradient method (abbreviated to IA).We set the starting point x 0 = x δ , and the parameters to τ k = 1 for all k, γ = 0.5, and ϕ i = 2, i = 1, 2. For the impulse noise, we use a noise density of 10%.In the case where x k+1 i is not set to x δ i , we use the scalar root solver scipy.optimize.brenthon Python.Otherwise, the updates are in closed form.
See Fig. 5 for numerical results.By gradient norm, we mean dist(∂ F(x k ), 0).We observe that, as in the convex case, the Bregman Itoh-Abe method converges significantly faster than the Itoh-Abe method, as the former utilises the sparse structure of the problem.

Conclusion
In this paper, we propose to discretise the ISS flow with the Itoh-Abe discrete gradient.The resultant schemes exhibit a dissipative structure (3.6) related to the symmetrised Bregman distance of a function J .This generalises the discrete gradient method for gradient flows and can be viewed as a discrete gradient analogue to Bregman iterations.Building on previous studies of the Itoh-Abe discrete gradient method in the non-smooth, non-convex setting, we prove convergence guarantees of the Bregman Itoh-Abe discrete gradient method in such a setting.
We consider numerical examples motivated by linear systems and searching for sparse solutions, as well as a nonconvex image denoising example.These results indicate that for sparse reconstructions, popular iterative solvers such as the SOR method can be significantly sped up by incorporating a Bregman step.
Future work is dedicated to proving convergence rates for the Bregman Itoh-Abe methods, and to compare the scheme to related methods such as the sparse Kaczmarz method [33].Furthermore, we will study corresponding inverse scale space schemes using other discrete gradients, such as the mean value discrete gradient.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.unmodified Bregman Itoh-Abe method (3.7) fails to converge to a limit set of stationary points.
Example A.1 As a starting point, we consider a C 2 -smooth objective function W : R 2 → R as described by Curry [17], for which the trajectory of a Euclidean gradient flow spirals along a gully, asymptotically trending towards the unit circle Denote by [ xk , ỹk ] T the iterates from the (standard) Itoh-Abe method for W with (arbitrary) time steps τ ≡ 1 and starting point [ x0 , ỹ0 ] T in the gully.As each update of the method consists of moving to a local point of descent, the iterates of this method will also remain in this gully and spiral towards S 1 .
We assume that there exists compact, disjoint sets A, B ⊂ R 2 with nonempty interior, such that B ∩ S 1 = ∅, and such that most of the iterates ([ xk , ỹk ] T ) k∈N are in A. That is, for all k ∈ N, the set {i ≤ k : [x i , ỹi ] T ∈ A} is strictly greater than the set {i ≤ k : [x i , ỹi ] T / ∈ A}.Such a set A exists as long as the iterates are spiralling around the unit circle at approximately a steadily decreasing rate.
Next we define an objective function V : R 3 → R as follows.Let f A , f B : R 2 → R be smooth support functions of A and B, respectively, with disjoint support sets.That is, f A (x, y) ∈ [0, 1], f B (x, y) ∈ [0, 1], and

Appendix B 1 -regularised sparse SOR method
In what follows, we describe the update rule for the Bregman Itoh-Abe method with V given in (5.5) and J given in (5.3).
Denote by xk+1 i the standard SOR update (5.2) for the ith coordinate.Then the 1 -regularised sparse SOR method can be expressed as follows.

Definition 2 . 3 (Example 2 . 4
Bregman distance) Given p ∈ ∂ J (x), the Bregman distance between x and y is given byD p J (x, y) = J (y) − J (x) − p, y − x .We refer to J as the corresponding Bregman function.If J (x) = x 2 /2, then D p J (x, y) = x−y 2 /2, i.e. the square of the Euclidean distance.

Proposition 2 . 14 (
Stationary versus minimal) Let V : R n → R be a convex function and C ⊂ R n be closed and convex.Then x * ∈ R n is a Clarke stationary point restricted to C if and only if x * ∈ arg min x∈C V (x).Proof This follows directly from [31, Theorem 4.19], noting that for convex functions the Clarke directional derivative coincides with the classical directional derivative [16, Proposition 2.2.7].

Assumption 3 . 1
(a) The function V : R n → R is locally Lipschitz continuous and bounded below.(b) x * ∈ R n is a Clarke stationary point of V restricted to

Remark 4 . 3
If the update is stationary, i.e. x k+1 i = x k i , then the subgradient update p k+1 i is unique only up to the choice of subderivative

Lemma 4 . 4 (
Properties of scheme) Let V , J , and C satisfy Assumption 3.1, and let (x k , p k ) k∈N be iterates that solve (3.7) for time steps (τ k ) k∈N ⊂ [τ min , τ max ].Then the following properties hold.

Theorem 4 . 5 (
Stationarity guarantees) Let V , J , and C satisfy Assumption 3.1, and suppose the sequence of iterates (x k , p k ) k∈N solves (4.3) for time steps (τ k ) k∈N ⊂ [τ min , τ max ].Then all accumulation points x * ∈ S are Clarke stationary points restricted to C.

Fig. 1
Fig. 1 Comparison of SOR and sparse SOR methods, for Gaussian linear system without noise.Top: Convergence rate for relative objective, i.e. [V (x k ) − V * ]/[V (x 0 ) − V * ].Bottom: Support error with respect to iterates, i.e. proportion of indices i s.t.sgn(x k i ) = sgn(x * i )

Fig. 2
Fig.2Comparison of SOR and sparse SOR methods, for Gaussian linear system without noise, and binary ground truth.Top: Convergence rate for relative objective.Bottom: Support error with respect to iterates

Fig. 4
Fig.4 Comparison of SOR and sparse SOR methods, for 1 -regularised linear system with noise.Top: Convergence rate for relative objective.Bottom: Support error with respect to iterates

Fig. 5
Fig. 5 Comparison of BIA and IA methods, for student-t regularised image denoising.First: Convergence rate for relative objective.Second: Convergence rate for relative gradient norm.Third: Input data.Fourth: Reconstruction f A (x, y) • f B (x, y) = 0 for all [x, y] T ∈ R 2 ,and furthermore, f A (A) ≡ 1 and f B (B) ≡ 1.We then consider the optimisation problem min z≥0 V (x, y, z) := W (x, y) + z( f A (x, y) − f B (x, y)).
2015 and a Professor in 2018.Since 2011, she has been a fellow of Jesus College, University of Cambridge, where she is currently a Professor of applied mathematics with DAMTP, the Head of the Cambridge Image Analysis Group, the Director of the Cantab Capital Institute for Mathematics of Information, and the Director of the Engineering and Physical Sciences Research Council Center for Mathematical and Statistical Analysis of Multimodal Clinical Imaging.Her research interests include variational methods, partial differential equations, and machine learning for image analysis, image processing, and inverse imaging problems.