Bregman Itoh--Abe methods for sparse optimisation

TL;DR

This paper introduces Bregman Itoh--Abe methods, a class of discrete gradient algorithms for sparse regularisation in inverse problems, offering guaranteed convergence and improved efficiency by leveraging problem structure.

Contribution

It develops Bregman discrete gradient methods for inverse problems, providing convergence guarantees and demonstrating their effectiveness in sparse optimisation.

Findings

Methods achieve rapid convergence rates.

Algorithms are unconditionally dissipative.

Numerical results show effectiveness in convex and non-convex cases.

Abstract

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.

Equations

$$\underset{\vec{x}\in \mathcal{C}}{\min }V(\vec{x}),$$

$$\dot{\vec{p}}(t)=-\nabla V(\vec{x}(t)),\vec{p}(t)\in \partial J(\vec{x}(t)).$$

$$S=\left\{{\vec{x}}^{\ast }\in {\mathbb{R}}^n\text{mathchar}58\exists ({\vec{x}}^{k_j}{)}_{j\in \mathbb{N}}\text{mbox}s.t.{\vec{x}}^{k_j}\to {\vec{x}}^{\ast }\right\}.$$

$$\mathrm{sgn}(x)=\{\begin{array}{ll}1,& \text{mbox}ifx>0,\\ 0,& \text{mbox}ifx=0,\\ -1,& \text{mbox}ifx<0.\end{array}$$

$$\partial J(\vec{x})=\{\vec{p}\in {\mathbb{R}}^n\text{mathchar}58J(\vec{y})-J(\vec{x})-⟨\vec{p},\vec{y}-\vec{x}⟩\ge 0\text{mbox}forall\vec{y}\in {\mathbb{R}}^n\}.$$

$${\chi }_{\mathcal{C}}(\vec{x})\text{mathchar}58=\{\begin{array}{ll}0,& \text{mbox}if\vec{x}\in \mathcal{C},\\ +\infty ,& \text{mbox}else.\end{array}$$

$$D_J^p(\vec{x},\vec{y})=J(\vec{y})-J(\vec{x})-⟨\vec{p},\vec{y}-\vec{x}⟩.$$

$$D_J^{\mathrm{symm}}(\vec{x},\vec{y})=D_J^q(\vec{y},\vec{x})+D_J^p(\vec{x},\vec{y})=⟨\vec{q}-\vec{p},\vec{y}-\vec{x}⟩.$$

$$D_J^p(\vec{x},\vec{y})\ge \frac{\mu }{2}\Vert \vec{y}-\vec{x}{\Vert }^2,D_J^{\mathrm{symm}}(\vec{x},\vec{y})\ge \mu \Vert \vec{y}-\vec{x}{\Vert }^2.$$

$$|V(\vec{y})-V(\vec{z})|\le L\Vert \vec{y}-\vec{z}\Vert .$$

$$V^o(\vec{x};\vec{d})=\underset{\vec{y}\to \vec{x},\lambda \downarrow 0}{\mathrm{limsup}}\frac{V(\vec{y}+\lambda \vec{d})-V(\vec{y})}{\lambda }.$$

$$\partial V(\vec{x})=\left\{\vec{p}\in {\mathbb{R}}^n\text{mathchar}58V^o(\vec{x};\vec{d})\ge ⟨\vec{d},\vec{p}⟩\text{mbox}forall\vec{d}\in {\mathbb{R}}^n\right\}.$$

$$\partial V(\vec{y})\subset \partial V(\vec{x})+B_{\epsilon }(\vec{0}),\text{mbox}forall\vec{y}\in B_{\delta }(\vec{x}).$$

$$\dot{\vec{x}}(t)=-({\nabla }^{2}J(\vec{x}(t)){)}^{-1}\nabla V(\vec{x}(t)),$$

$$\frac{d}{dt}V(\vec{x}(t))$$

$$=-⟨\dot{\vec{x}}(t),{\nabla }^{2}J(\vec{x}(t))\dot{\vec{x}}(t)⟩\le -\mu \Vert \dot{\vec{x}}(t){\Vert }^2.$$

$${\vec{p}}^{k+1}={\vec{p}}^k-{\tau }_k\nabla V({\vec{x}}^{k+1}),{\vec{p}}^{k+1}\in \partial J({\vec{x}}^{k+1})$$

$${\vec{x}}^{k+1}=\underset{\vec{x}\in {\mathbb{R}}^n}{\mathrm{argmin}}V(\vec{x})+\frac{1}{{\tau }_k}{D}_J^{{\vec{p}}^k}({\vec{x}}^k,\vec{x}).$$

$$V({\vec{x}}^{k+1})-V({\vec{x}}^k)$$

$$\le -\frac{\mu }{2{\tau }_k}\Vert {\vec{x}}^k-{\vec{x}}^{k+1}{\Vert }^2.$$

$${\vec{p}}^{k+1}={\vec{p}}^k-{\tau }_k\nabla V({\vec{x}}^k),{\vec{p}}^{k+1}\in \partial J({\vec{x}}^{k+1})$$

$${\vec{x}}^{k+1}=\underset{\vec{x}\in {\mathbb{R}}^n}{\mathrm{argmin}}V({\vec{x}}^k)+⟨\nabla V({\vec{x}}^k),\vec{x}-{\vec{x}}^k⟩+\frac{1}{{\tau }_k}{D}_J^{{\vec{p}}^k}({\vec{x}}^k,\vec{x}).$$

$$V(\vec{x})=\frac{1}{2}\Vert A\vec{x}-\vec{b}{\Vert }^2.$$

$$\underset{\vec{x}\in {\mathbb{R}}^n}{\min }\{J(\vec{x})\text{mbox}s.t.A\vec{x}=\vec{b}\}.$$

$$⟨\overline{\nabla }V(\vec{x},\vec{y}),\vec{y}-\vec{x}⟩$$

$$\underset{\vec{y}\to \vec{x}}{\lim }\overline{\nabla }V(\vec{x},\vec{y})$$

$$\overline{\nabla }V(\vec{x},\vec{y})=\left(\begin{array}{c}\frac{V\left(y_1,x_2,\dots ,x_n\right)-V(\vec{x})}{y_1-x_1}\\ \frac{V\left(y_1,y_2,x_3,\dots ,x_n\right)-V\left(y_1,x_2,\dots ,x_n\right)}{y_2-x_2}\\ ⋮\\ \frac{V(\vec{y})-V\left(y_1,\dots ,y_{n-1},x_n\right)}{y_n-x_n}\end{array}\right),$$

$${\vec{p}}^{k+1}$$

$$V({\vec{x}}^k)-V({\vec{x}}^{k+1})$$

$$=\frac{1}{{\tau }_k}⟨{\vec{x}}^k-{\vec{x}}^{k+1},{\vec{p}}^k-{\vec{p}}^{k+1}⟩$$

$$=\frac{1}{{\tau }_k}{D}_J^{\mathrm{symm}}({\vec{x}}^k,{\vec{x}}^{k+1})$$

$$\ge \frac{\mu }{{\tau }_k}\Vert {\vec{x}}^k-{\vec{x}}^{k+1}{\Vert }^2.$$

$$\partial J(\vec{x})=\prod _{i=1}^n\partial {\chi }_{[l_i,u_i]}(x_i)+\partial j_i(x_i).$$

Full text

Keywords.   Non-convex optimisation, non-smooth optimisation, Bregman iteration, inverse scale space, geometric numerical integration, discrete gradient methods

AMS subject classifications.   49M37, 49Q15, 65K10, 90C26

††journal: Journal of Mathematical Imaging and Vision

1 Introduction

We consider the constrained optimisation problem

[TABLE]

for an objective function $V\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\mathbb{R}$ and constraint $\mathcal{C}\subset\mathbb{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\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\overline{\mathbb{R}}$. The ISS flow is given by

[TABLE]

Here, $\partial J$ is the convex subdifferential of $J$, to be defined in Section 2, and $\overline{\mathbb{R}}\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{R}\cup\{\pm\infty\}$. The term inverse scale space flow goes back to Scherzer & Groetsch (sch01, ). This flow is typically derived as the continous-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 Section 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 (rii18smooth, ; gri17, ; rii18nonsmooth, ; rin18, ), 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 (rii18nonsmooth, ). For smooth problems, the theoretical convergence rates of the discrete gradient methods match those of explicit gradient descent and coordinate descent (rii18smooth, ). 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. (miy18, ) 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.

1.1 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. (osh05, ) in 2005 for total variation-based image denoising, representing an extension of the Bregman proximal algorithm (cen92, ; eck93, ; kiw97, ; teb92, ) to non-smooth Bregman functions. Subsequently the ISS flow was derived and analysed by Burger et al. (bur06, ; bur07, ; bur13iss, ; bur16bregman, ). Since then, researchers have studied the ISS flow with applications to generalised spectral analysis in a nonlinear setting, i.e. by Burger et al. (bur16, ), Gilboa et al. (gil16, ), and Schmidt et al. (sch18, ). The Bregman method has been studied for $\ell^{1}$-regularisation and compressed sensing by Goldstein & Osher (gol09, ) and Yin et al. (yin08, ), and extended to primal-dual algorithms by Zhang et al. (zha11, ).

The linearised Bregman method was proposed by Yin et al. (yin08, ) for applications to $\ell^{1}$-regularisation and compressed sensing, and further studied in this setting by Cai et al. (cai09, ), and Dong et al. (don10, ). An extension for non-convex problems was proposed by Benning et al. (ben17, ), proving global convergence for functions that satisfy the Kurdyka–Łojasiewicz property. Lorenz et al. (lor14, ; sch18kaczmarz, ) 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 genereal convergence results by Lorenz et al. (lor14sfp, ). For further details on Bregman iterations and linearised Bregman methods, we refer to (ben18, ).

We review papers on discrete gradient methods for optimisation based on gradient flows. Grimm et al. (gri17, ) 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. (rii18smooth, ) 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. (rin18, ) 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. (rii18nonsmooth, ) 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. (cel18, ) extended the Itoh–Abe discrete gradient method to optimisation problems defined on Riemannian manifolds. Hernández-Solano et al. (her15, ) 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. (su16, ) and Wibisono et al. (wib16, ) study second-order ODEs based on the continuous-time limit of Nesterov’s accelerated gradient descent (nes83, ), in order to gain insight into the acceleration phenomenon. In this setting, Wilson et al. (wil16, ) provided a framework for Lyapunov analysis of optimisation schemes and continuous-time dynamics, and Betancourt et al. (bet18, ) proposed a framework for symplectic integration for optimisation. On a related note, Scieur et al. (sci17, ) 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. (mad18, ), with the aim of ensuring linear convergence for a larger group of functions than classical methods. Other numerical integration methods, such as implicit Runge-Kutta methods, where energy dissipation is ensured under mild time step restrictions (hai13, ), and explicit stabilised methods for solving strongly convex problems (eft18, ). Another example is the study of gradient flows in metric spaces and minimising movement schemes (amb08, ), concerning gradient flow trajectories under other measures of distance, such as the Wasserstein metric (san17, ).

1.2 Structure and contributions

The paper is structured as follows. For the remainder of Section 1, we define mathematical notation. In Section 2, we provide preliminary material for convex and non-convex analysis, and introduce Bregman distances. In Section 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 Section 4, we discuss particular examples of Bregman discrete gradient methods. In Section 5, we present results from numerical experiments. In Section 6, we conclude.

1.3 Notation and preliminaries

Throughout the paper, we denote by $\|\cdot\|$ the $\ell^{2}$-norm, i.e. $\|\vec{x}\|^{2}=\sum_{i=1}^{n}|x_{i}|^{2}$. For $p\in[1,+\infty]$, $\|\cdot\|_{p}$ denotes the usual $\ell^{p}$-norm. For $\varepsilon>0$ and $\vec{x}\in\mathbb{R}^{n}$, we denote by $B_{\varepsilon}(\vec{x})$ the open ball $\{\vec{y}\in\mathbb{R}^{n}\;\mathrel{\mathop{\mathchar 58\relax}}\;\|\vec{y}-\vec{x}\|<\varepsilon\}$. We denote by $\{\vec{e}^{1},\ldots,\vec{e}^{n}\}$ the standard coordinate vectors in $\mathbb{R}^{n}$. We define the extended reals as $\overline{\mathbb{R}}\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{R}\cup\{\pm\infty\}$.

For a sequence $(\vec{x}^{k})_{k\in\mathbb{N}}\subset\mathbb{R}^{n}$, we denote by $S$ its limit set, i.e. the set of accumulation points,

[TABLE]

For a matrix $A\in\mathbb{R}^{m\times n}$, we denote by $\vec{a}^{i}\in\mathbb{R}^{n}$ its $i$th row, and $A^{T}$ its transpose. Accordingly, $a^{i}_{j}$ refers to the $j$th element of the $i$th column of $A$.

For $x\in\mathbb{R}$, the sign operator $\operatorname{sgn}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}\to\{\pm 1,0\}$ is defined as

[TABLE]

2 Preliminaries for convex and non-convex analysis

In this section, we provide preliminary material on nonsmooth analysis. In Section 2.1, we cover convex analysis and Bregman distances, while in Section 2.2 we cover nonconvex, nonsmooth analysis in the Clarke subdifferential framework.

2.1 Convex analysis

We consider functions $J\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\overline{\mathbb{R}}$ that are convex, proper, and lower-semicontinuous (see (eke99, ; roc15, ) for details on this class of functions).

Definition 2.1** (Effective domain)**

The effective domain of a function $J\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\overline{\mathbb{R}}$ is defined as $\operatorname{dom}(J)=\{\vec{x}\in\mathbb{R}^{n}\;\mathrel{\mathop{\mathchar 58\relax}}\;J(\vec{x})<\infty\}$.

Definition 2.2** (Subgradients and subdifferentials)**

The subdifferential of a convex function $J\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\overline{\mathbb{R}}$ at $\vec{x}\in\mathbb{R}^{n}$ is the set of vectors

[TABLE]

Vectors in $\partial J(\vec{x})$ are called subgradients of $J$ at $\vec{x}$.

For $i=1,\ldots,n$, we denote by $[\partial J(\vec{x})]_{i}$ the projection of $\partial J(\vec{x})$ onto the $i$th coordinate, i.e. $[\partial J(\vec{x})]_{i}=\{p_{i}\;\mathrel{\mathop{\mathchar 58\relax}}\;\vec{p}\in\partial J(\vec{x})\}$.

For the theory of convex functions and their subdifferentials, see (roc15, ).

We consider the characteristic function $\chi_{\mathcal{C}}$ of a convex, closed set $\mathcal{C}\subset\mathbb{R}^{n}$, defined by

[TABLE]

This is a convex, proper, lower-semicontinuous function, and $\vec{0}\in\partial\chi_{\mathcal{C}}(\vec{x})$ for all $\vec{x}\in\operatorname{dom}(\chi_{\mathcal{C}})=\mathcal{C}$.

We can now define the generalised Bregman distance (bre67, ) of a convex function $J$.

Definition 2.3** (Bregman distance)**

Given $\vec{p}\in\partial J(\vec{x})$, the Bregman distance between $\vec{x}$ and $\vec{y}$ is given by

[TABLE]

We refer to $J$ as the corresponding Bregman function.

Example 2.4

If $J(\vec{x})=\|\vec{x}\|^{2}/2$, then $D^{p}_{J}(\vec{x},\vec{y})=\|\vec{x}-\vec{y}\|^{2}/2$, i.e. the square of the Euclidean distance.

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 $\vec{p}\in\partial J(\vec{x})$ and $\vec{q}\in\partial J(\vec{y})$, the symmetric Bregman distance between $\vec{x}$ and $\vec{y}$ is given by

[TABLE]

Definition 2.6** ($\mu$-convexity)**

A proper, lower-semicontinuous convex function $J\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\overline{\mathbb{R}}$ is $\mu$-convex for $\mu\geq 0$ if either of the following (equivalent) conditions hold.

1. (i)

The function $J(\cdot)-\dfrac{\mu}{2}\|\cdot\|^{2}$ is convex. 2. (ii)

$J(\vec{y})-J(\vec{x})-\langle\vec{p},\vec{y}-\vec{x}\rangle\geq\tfrac{\mu}{2}\|\vec{y}-\vec{x}\|^{2}$* for all $\vec{x},\vec{y}\in\mathbb{R}^{n}$, $\vec{p}\in\partial J(\vec{x})$.*

If $\mu>0$, $J$ is said to be strongly convex.

For strongly convex Bregman functions, the following property of Bregman distances is immediate.

Proposition 2.7** (Bregman distance under $\mu$-convexity)**

If $J$ is $\mu$-convex, then for all $\vec{x},\vec{y}\in\mathbb{R}^{n}$,

[TABLE]

2.2 Non-convex subdifferential analysis

We summarise the main concepts of the Clarke subdifferential for locally Lipschitz continuous, non-smooth, non-convex functions $V\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\mathbb{R}$, and refer to (cla90, ) for further details.

Definition 2.8** (Lipschitz continuity)**

$V$* is Lipschitz of rank $L$ near $\vec{x}$ if there exists $\varepsilon>0$ such that for all $\vec{y},\vec{z}\in B_{\varepsilon}(\vec{x})$, one has*

[TABLE]

$V$* is locally Lipschitz continuous if the above property holds for all $\vec{x}\in\mathbb{R}^{n}$.*

Definition 2.9** (Clarke directional derivative)**

For $V$ Lipschitz near $\vec{x}$ and for a vector $\vec{d}\in\mathbb{R}^{n}$, the Clarke directional derivative is given by

[TABLE]

Definition 2.10** (Clarke subdifferential)**

Let $V$ be locally Lipschitz continuous and $\vec{x}\in\mathbb{R}^{n}$. The Clarke subdifferential of $V$ at $\vec{x}$ is given by

[TABLE]

An element of $\partial V(\vec{x})$ is called a Clarke subgradient.

The Clarke subdifferential was introduced by Clarke in (cla73, ). It is well-defined for locally Lipschitz functions, coincides with the standard subdifferential for convex functions (cla90, , Proposition 2.2.7), and coincides with the derivative at points of strict differentiability (cla90, , Proposition 2.2.4). We additionally state two useful results, both of which can be found in Chapter 2 of (cla90, ).

Proposition 2.11** (Properties of Clarke subdifferential)**

Suppose $V$ is locally Lipschitz continuous. Then

1. (i)

$\partial V(\vec{x})$* is nonempty, convex and compact, and if $V$ is Lipschitz of rank $L$ near $\vec{x}$, then $\partial V(\vec{x})\subseteq B_{L}(\vec{0})$.* 2. (ii)

$\partial V(\vec{x})$* is outer semicontinuous. That is, for all $\varepsilon>0$, there exists $\delta>0$ such that*

[TABLE]

In this paper, we will consider constrained optimisation problems $\operatorname*{arg\,min}_{\vec{x}\in\mathcal{C}}V(\vec{x})$ for a closed, convex subset $\mathcal{C}$. We assume in the convergence analysis that $\mathcal{C}$ is box constraints of the form $\mathcal{C}=\otimes_{i=1}^{n}[l_{i},u_{l}]$. To have a formal notion of first-order optimality in this setting, we define the tangent cone.

Definition 2.12** (Tangent cone)**

For a convex, closed set $\mathcal{C}\subset\mathbb{R}^{n}$, the tangent cone at $\vec{x}\in\mathcal{C}$, denoted by $T(\mathcal{C},\vec{x})$, is the closure of the set of vectors $\vec{v}\in\mathbb{R}^{n}$ such that there is $\delta>0$ such that for all $\varepsilon\in(0,\delta)$, $\vec{x}+\varepsilon\vec{v}\in\mathcal{C}$.

Definition 2.13** (Clarke stationarity)**

Let $V\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\mathbb{R}$ be a locally Lipschitz continuous function and $\mathcal{C}\subset\mathbb{R}^{n}$ a closed, convex set. A point $\vec{x}^{*}\in\mathbb{R}^{n}$ is a Clarke stationary point of $V$ restricted to $\mathcal{C}$ if $V^{o}(\vec{x}^{*};\vec{d})\geq 0$ for all $\vec{d}\in T(\mathcal{C},\vec{x}^{*})$.

As the following proposition shows, the above definition generalises minimisers for convex functions, and Clarke stationary points in the unrestricted case.

Proposition 2.14** (Stationary versus minimal)**

Let $V\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\mathbb{R}$ be a convex function and $\mathcal{C}\subset\mathbb{R}^{n}$ be closed and convex. Then $\vec{x}^{*}\in\mathbb{R}^{n}$ is a Clarke stationary point restricted to $\mathcal{C}$ if and only if $\vec{x}^{*}\in\operatorname*{arg\,min}_{\vec{x}\in\mathcal{C}}V(\vec{x})$.

Proof

This follows directly from (jah07, , Theorem 4.19), noting that for convex functions the Clarke directional derivative coincides with the classical directional derivative (cla90, , Proposition 2.2.7).

3 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.

3.1 Inverse scale space flow and Bregman methods

For a convex function $J\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\overline{\mathbb{R}}$, objective function $V\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\mathbb{R}$ and starting points $\vec{x}(0)=\vec{x}^{0}\in\mathbb{R}^{n}$, $\vec{p}(0)\in\partial J(\vec{x}^{0})$, the ISS flow is the dissipative differential system given by (1.2). If $J$ were twice continuously differentiable and $\mu$-convex, then (1.2) could be rewritten as

[TABLE]

and the energy $V(\vec{x}(t))$ would dissipate over time as

[TABLE]

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

[TABLE]

which can be rewritten as

[TABLE]

From (3.1), we see that the Bregman method is dissipative, as

[TABLE]

Similarly, the linearised Bregman method is derived by forward Euler discretisation of (1.2), and is given by

[TABLE]

or equivalently

[TABLE]

The ISS flow and Bregman methods are considered for solving ill-conditioned linear systems $A\vec{x}=\vec{b}$, with the objective function

[TABLE]

In this setting, iterates of both the Bregman method and the linearised Bregman method converge (ben18, ; lor14sfp, ) to a solution of

[TABLE]

Furthermore, applications of the ISS flow include image denoising with reduced contrast-loss and staircasing effects (osh05, ), recovering eigenfunctions (sch18, ), and identifying sparse or low-rank structures (yin08, ).

We make the following assumptions on the the objective function $V$, the constraints $\mathcal{C}$, as well as the Bregman function $J$.

Assumption 3.1

1. a)

The function $V\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\mathbb{R}$ is locally Lipschitz continuous and bounded below. 2. b)

$\vec{x}^{*}\in\mathbb{R}^{n}$* is a Clarke stationary point of $V$ restricted to $\mathcal{C}$ if and only if for all coordinate vectors $\vec{e}^{i}$, we have $V^{o}(\vec{x}^{*};\vec{e}^{i}),V^{o}(\vec{x}^{*};-\vec{e}^{i})\geq 0$.* 3. c)

The set $\mathcal{C}\subset\mathbb{R}^{n}$ consists of coordinate-wise box constraints, i.e. $\mathcal{C}=\otimes_{i=1}^{n}[l_{i},u_{i}]$. 4. d)

The function $J\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\overline{\mathbb{R}}$ is proper, lower-semicontinuous, and $\mu$-convex with $\mu>0$. Furthermore, $J(\vec{x})=\sum_{i=1}^{n}j_{i}(x_{i})+\chi_{[l_{i},u_{i}]}(x_{i})$, where $j_{i}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}\to\overline{\mathbb{R}}$ are convex and Lipschitz continuous on $[l_{i},u_{i}]$ for each $i$.

3.2 The Bregman discrete gradient method

In what follows, we define discrete gradients, and propose the Bregman discrete gradient method by discretising the ISS flow.

Definition 3.2** (Discrete gradient)**

Let $V$ be a continuously differentiable function. A discrete gradient is a continuous map $\overline{\nabla}V\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ such that for all $\vec{x},\vec{y}\in\mathbb{R}^{n}$,

[TABLE]

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 (hai06, ; mcl01, ) for an introduction. Discrete gradients are tools for solving first-order ODEs that preserve energy conservation laws, dissipation laws, and Lyapunov functions (gon96, ; ito88, ; mcl99, ; qui96, ).

The Itoh–Abe discrete gradient (ito88, ) (also known as the coordinate increment discrete gradient) is given by

[TABLE]

where $0/0$ is interpreted as $\partial_{i}V(\vec{x})$.

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. (rin18, )), or unavailable (e.g. (rii18nonsmooth, )), Itoh–Abe discrete gradient methods are useful for ensuring convergence guarantees and competitive convergence rates (rii18smooth, ) without requiring derivatives.

We propose the Bregman discrete gradient method as follows. For a starting point $\vec{x}^{0}\in\mathbb{R}^{n}$, subgradient $\vec{p}^{0}\in\partial J(\vec{x}^{0})$, and time steps $(\tau_{k})_{i\in\mathbb{N}}>0$, solve for $k=0,1,\ldots$,

[TABLE]

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).

[TABLE]

Furthermore, if we plug in $J(\vec{x})=\|\vec{x}\|^{2}/2$, we recover the discrete gradient method for gradient flows (rii18smooth, ). Observe that the dissipative structure holds for arbitrary positive time steps $\tau_{k}>0$. We thus only require that $\tau_{k}\in[\tau_{\min},\tau_{\max}]$ for arbitrary bounds $\tau_{\min},\tau_{\max}>0$.

By Assumption 3.1 d), the subdifferential of $J$ is separable in the coordinates, i.e.

[TABLE]

It follows that solving the Bregman Itoh–Abe equation (3.5) corresponds to successively solving $n$ scalar inclusions, i.e. given $\vec{x}^{k}$, and $\vec{p}^{k}$, we want $\vec{x}^{k+1}$ and $\vec{p}^{k+1}$ that solve

[TABLE]

Here $\vec{y}^{k,i}$ denotes $[x^{k+1}_{1},\ldots,x^{k+1}_{i},x^{k}_{i+1},\ldots,x^{k}_{n}]^{T}$. For a choice of $v^{k}_{i}\in[\partial V(\vec{y}^{k,i-1})]_{i}$, if $p^{k}_{i}-\tau_{k,i}v^{k}_{i}\in[\partial J(\vec{y}^{k,i-1})]_{i}$, then we consider $x^{k+1}_{i}=x^{k}_{i}$ and $p^{k+1}_{i}=p^{k}_{i}-\tau_{k,i}v^{k}_{i}$ an admissible update.

Remark 3.3

While discrete gradients are not defined for non-smooth 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 (cla90, ). Furthermore, note that the properties used to derive the dissipative structure (3.6) holds in this setting.

4 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 $(\vec{x}^{k})_{k\in\mathbb{N}}$ from the Bregman Itoh–Abe method are Clarke stationary.

4.1 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 adress these questions. The corresponding analysis generalises that from the previous papers on discrete gradient methods for optimisation (rii18smooth, ; rii18nonsmooth, ).

We first show that an update (not necessarily unique) always exists.

Lemma 4.1** (Existence)**

For any $\tau>0$, $\vec{x}^{k}\in\mathbb{R}^{n}$, and $\vec{p}^{k}\in\partial J(\vec{x}^{k})$, there exists an update $(\vec{x}^{k+1},\vec{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\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}\to\mathbb{R}$, $j\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}\to\overline{\mathbb{R}}$. For $x\in\mathbb{R}$ and $p\in\partial j(x)$ we either want $y\neq x$ such that

[TABLE]

or $y=x$ and $p-\tau w\in\partial j(x)$, for some $w\in\partial 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-\tau v^{o}(x;1)>p$ and $p\in\partial j(x)$, we deduce that $p-\tau v^{o}(x;1)>\partial j(x)$. By the outer semicontinuity of subdifferentials and definition of Clarke directional derivatives, there is $\delta>0$ such that

[TABLE]

On the other hand, as $v$ is bounded below,

[TABLE]

is bounded below for all $y\in[x+\delta,+\infty)$, while by $\mu$-convexity of $j$, we have $\partial j(y)\geq\partial j(x)+\mu(y-x)$ for all $y\in[x+\delta,+\infty)$. Hence, there is $r\gg 0$ such that

[TABLE]

By continuity of $v$, and by outer semicontinuity of subdifferentials, it follows that there exists $y\in(x+\delta,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 (rii18smooth, , Section 5).

Proposition 4.2** (Uniqueness)**

For $\tau>0$, $\vec{x}^{k}\in\mathbb{R}^{n}$, and $\vec{p}^{k}\in\partial J(\vec{x}^{k})$, the update $(\vec{x}^{k+1},\vec{p}^{k+1})$ to (3.7) is unique in the following cases.

1. (i)

$V$* is convex.* 2. (ii)

$V(\cdot)+\eta\|\cdot\|^{2}$* is convex for $\eta>0$, and $\tau<\mu/\eta$.*

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 $\vec{y}^{k,i}$ must satisfy

[TABLE]

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(\cdot)+\eta\|\cdot\|^{2}$ is convex and $\tau<\mu/\eta$. We rewrite (4.2) as

[TABLE]

Since $\tau\eta<\mu$, the function $J(\cdot)-\tau\eta\|\cdot\|^{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 $V(\cdot)+\eta\|\cdot\|^{2}$ is convex, the difference quotient term $\tfrac{V(\vec{y}^{k,i})-V(\vec{y}^{k,i-1})}{x^{k+1}_{i}-x^{k}_{i}}+\eta x^{k+1}_{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.

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 $v_{i}\in[\partial V(\vec{y}^{k,i-1})]_{i}$.

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 (roc85, ) and prox-regularity (pol96, ). 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 $(\vec{x}^{k},\vec{p}^{k})_{k\in\mathbb{N}}$, and generalises Lemmas 3.3, 3.4 in (rii18nonsmooth, ).

Lemma 4.4** (Properties of scheme)**

Let $V$, $J$, and $\mathcal{C}$ satisfy Assumption 3.1, and let $(\vec{x}^{k},\vec{p}^{k})_{k\in\mathbb{N}}$ be iterates that solve (3.7) for time steps $(\tau_{k})_{k\in\mathbb{N}}\subset[\tau_{\min},\tau_{\max}]$. Then the following properties hold.

1. (i)

$V(\vec{x}^{k+1})\leq V(\vec{x}^{k})$. 2. (ii)

$\lim_{k\to\infty}\|\vec{x}^{k+1}-\vec{x}^{k}\|=0$. 3. (iii)

If $V$ is coercive, then there exists a convergent subsequence of $(\vec{x}^{k},\vec{p}^{k})_{k\in\mathbb{N}}$. 4. (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(\vec{x}^{k}))_{k\in\mathbb{N}}$ is decreasing, $V(\vec{x}^{k})\to V^{*}$ for some limit $V^{*}$. Therefore, by (3.6),

[TABLE]

This implies property (ii).

Properties (iii) and (iv) follow from (i) and (ii) and are proven respectively in (rii18nonsmooth, , Lemma 3.3 and Lemma 3.4).

4.2 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 $\vec{x}^{0}$, $\vec{p}^{0}\in\partial J(\vec{x}^{0})$, and $k\in\mathbb{N}$, update

[TABLE]

Observe that since $\vec{0}\in\partial J(\vec{x})$ for all $\vec{x}\in\mathcal{C}$, we still have $\vec{p}^{k+1}\in\partial J(\vec{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, $\partial\chi_{[l_{i},u_{i}]}$ is unbounded). Otherwise, the iterates might get stuck at a nonstationary point for several iterations, since $x^{k}_{i}$ is unable to leave $\{l_{i},u_{i}\}$ until the subgradient update in $\partial\chi_{[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.

4.3 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 nonstationary. We note that this also holds for the original Bregman Itoh–Abe method if the iterates $(\vec{x}^{k})_{k\in\mathbb{N}}$ converge to a unique limit.

Theorem 4.5** (Stationarity guarantees)**

Let $V$, $J$, and $\mathcal{C}$ satisfy Assumption 3.1, and suppose the sequence of iterates $(\vec{x}^{k},\vec{p}^{k})_{k\in\mathbb{N}}$ solves (4.3) for time steps $(\tau_{k})_{k\in\mathbb{N}}\subset[\tau_{\min},\tau_{\max}]$. Then all accumulation points $\vec{x}^{*}\in S$ are Clarke stationary points restricted to $\mathcal{C}$.

Proof

Let $\vec{x}^{*}\in S$ and consider a convergent subsequence $(\vec{x}^{k_{j}})_{j\in\mathbb{N}}$. We want to show for each basis vector $\vec{e}^{i}$ that either $V^{o}(\vec{x}^{*};\vec{e}^{i})\geq 0$ or $x^{*}_{i}=u_{i}$, and analogously that either $V^{o}(\vec{x}^{*};-\vec{e}^{i})\geq 0$ or $x^{*}_{i}=l_{i}$. As the arguments are equivalent, we only prove the first case.

Suppose for contradiction that $V^{o}(\vec{x}^{*};\vec{e}^{i})<-\eta$ for some $\eta>0$, and that $x^{*}_{i}<u_{i}$. By the definition of the Clarke directional derivative, Definition 2.9, there are $\varepsilon,\delta>0$ such that for all $\vec{x}\in B_{\varepsilon}(\vec{x}^{*})$ and $\lambda\in(0,\delta)$, we have

[TABLE]

Since $\vec{x}^{k_{j}}\to\vec{x}^{*}$ and $\|\vec{x}^{k_{j}+1}-\vec{x}^{k_{j}}\|\to 0$, for each $N\in\mathbb{N}$ there exists $K$ such that for all $j\geq K$, we have $\vec{x}^{k}\in B_{\varepsilon}(\vec{x}^{*})$ and $\|\vec{x}^{k}-\vec{x}^{k+1}\|<\delta$ for $k=k_{j},k_{j}+1,\ldots,k_{j}+N$. By making $\varepsilon>0$ sufficiently small, we have $B_{\varepsilon}(x^{*}_{i})<u_{i}$. Furthermore, since $x^{k+1}_{i}\geq x^{k}_{i}$ for $k=k_{j},\ldots,k_{j}+N-1$, we deduce that the constraint component $\tilde{q}$ is zero. By (4.4), it follows that

[TABLE]

By Assumption 3.1, $\partial j_{i}$ is bounded on $U=B_{\varepsilon}(\vec{x}^{*})\cap[l_{i},u_{i}]$. Since $p^{k,j}_{i},\ldots,p^{k_{j}+N}_{i}\in\partial j_{i}(U)$, we can choose $N$ such that $N\tau_{\min}\frac{\eta}{2}>\max\partial j_{i}(U)-\min\partial j_{i}(U)$ and arrive at a contradiction. Thus, $\vec{x}^{*}$ is a Clarke stationary point restricted to $\mathcal{C}$.

5 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

[TABLE]

where $A\in\mathbb{R}^{n\times 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 $\vec{\tau}^{k}$ coordinate-wise scaled by the diagonal of $A$, i.e. $\vec{\tau}^{k}=\tau/\mbox{diag}\,(A)=[\tau/a^{1}_{1},\ldots,\tau/a^{n}_{n}]$ for all $k\in\mathbb{N}$, and some $\tau>0$.

We first introduce some well-known coordinate descent schemes for solving linear systems, which Miyatake et al. (miy18, ) showed were equivalent to the Itoh–Abe discrete gradient method. The successive-over-relaxation (SOR) method (you71, ) updates each coordinate sequentially according to the rule

[TABLE]

where $\omega\in(0,2)$. For $\omega=1$, this is the Gauss-Seidel method (you71, ). The SOR method is equivalent to the Itoh–Abe discrete gradient method

[TABLE]

with $V$ given by (5.1) with the time steps $\tau_{i}=2\omega/\mathinner{\left((2-\omega)a^{i}_{i}\right)}$.

5.1 Sparse SOR method

We consider underdetermined linear systems and want to find sparse solutions $\vec{x}^{*}$. Hence we seek to apply the Bregman Itoh–Abe method (3.7) with objective function $V$ given by (5.1), and

[TABLE]

for $\gamma>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 $\tilde{x}^{k+1}_{i}$ the standard SOR coordinate update from $x^{k}_{i}$, (5.2). Furthermore, for $\vec{p}^{k}\in\partial J(\vec{x}^{k})$, we write $\vec{p}^{k}=\vec{x}^{k}+\gamma\vec{r}^{k}$, where $\vec{r}^{k}\in\partial\|\vec{x}^{k}\|_{1}$. Then $(x^{k+1}_{i},r^{k+1}_{i})$ are given in closed form as

[TABLE]

Here $S\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ denotes the shrinkage operator $S(\vec{x},\lambda)=\operatorname{sgn}(\vec{x})\max\{|\vec{x}|-\lambda,0\}$, applied elementwise to a vector $\vec{x}$.

5.2 Sparse, regularised SOR

If $\vec{b}=A\vec{x}^{\text{true}}+\vec{\delta}$, where $\vec{x}^{\text{true}}$ is the sparse ground truth and $\vec{\delta}$ is noise, then it may be necessary to regularise the objective function as well. Hence we consider the objective function

[TABLE]

for some regularisation parameter $\lambda>0$. The non-smoothness induced by $\|\cdot\|_{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(\vec{x})=\frac{1}{2}\|\vec{x}\|^{2}+\gamma\|\vec{x}\|_{1}$ and $J(\vec{x})=\frac{1}{2}\|\vec{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 include this in the appendix.

6 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 (miy18, ).

The explicit coordinate descent method (bec13, ; wri15, ) for minimising $V$ is given by

[TABLE]

where $\alpha_{i}>0$ is the time step. As mentioned in (wri15, ), the SOR method is also equivalent to the coordinate descent method with $V$ given by (5.1) and the time step $\alpha_{i}=\omega/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 Section 5.1. We may compare this to a Bregman linearised coordinate descent scheme,

[TABLE]

One can verify that this scheme is equivalent to (6.1) for the parameters

[TABLE]

7 Numerical examples

In this section, we present numerical results for the schemes described in Section 5.

7.1 Sparse SOR

We construct a matrix $A\in\mathbb{R}^{1024\times 1024}$ from independent standard (zero mean, unit variance) Gaussian draws, and construct the sparse ground truth $\vec{x}^{\text{true}}$ by choosing $10\%$ of the indices at random determined by uniform draws on the unit interval. We then solve the problem

[TABLE]

where $\vec{b}=A\vec{x}^{\text{true}}$. We compare the SOR method ($J(\vec{x})=\|\vec{x}\|^{2}/2$) and the BSOR method ($J(\vec{x})=\|\vec{x}\|^{2}/2+\gamma\|\vec{x}\|_{1}$), where $\gamma=1$. We set time steps to $\vec{\tau}=2/\mbox{diag}\,(A)$, corresponding to the Gauss-Seidel method. See Figure 7.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 Figure 7.2.

7.2 Sparse, regularised SOR

We construct $A\in\mathbb{R}^{1024\times 1024}$ and $\vec{x}^{\text{true}}$ as in the previous subsection. However, we add noise to the data, i.e. $\vec{b}=A\vec{x}^{\text{true}}+\vec{\delta}$, where $\vec{\delta}$ is independent Gaussian noise with a standard deviation of $0.1\|A\vec{x}^{\text{true}}\|_{\infty}$. Since the added noise destroys the sparsity structure of $A^{-1}\vec{b}$, the sparse SOR method fails to improve the convergence rate. The results for $V(\vec{x})=\|A\vec{x}-\vec{b}\|^{2}/2$ are in Figure 7.3.

We therefore include regularisation in the objective function of the form

[TABLE]

where $\lambda=100$, and with initialisation $\vec{x}^{0}$ constructed by random, independent Gaussian draws. The results are visualised in Figure 7.4.

7.3 Student-t regularised image denoising

We consider a nonconvex image denoising model, previously presented in (och14, ), given by

[TABLE]

Here $\{K_{i}\}_{i=1}^{N}$ is a collection of linear filters, $(\varphi_{i})_{i=1}^{N}\subset[0,\infty)$ are coefficients, $\Phi\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\mathbb{R}$ is the nonconvex function based on the t-student distribution, defined as

[TABLE]

and $x^{\delta}$ is an image corrupted by impulse noise (salt & pepper noise).

As impulse noise only affects a small subset of pixels, we use the data fidelity term $\vec{x}\mapsto\|\vec{x}-\vec{x}^{\delta}\|_{1}$ to promote sparsity of $\vec{x}^{*}-\vec{x}^{\delta}$ for $\vec{x}^{*}\in\operatorname*{arg\,min}F(\vec{x})$. As linear filters, we consider the simple case of finite difference approximations to first-order derivatives of $\vec{x}$. We note that by applying a gradient flow to this regularisation function, we observe a similarity to Perona-Malik diffusion (per90, ).

We consider a Bregman Itoh–Abe method (abbreviated to BIA) for solving $\min_{\vec{x}\in\mathbb{R}^{n}}F(\vec{x})$ with the Bregman function

[TABLE]

to account for the sparsity of the residual $\vec{x}^{*}-\vec{x}^{\delta}$, and compare the method to the regular Itoh–Abe discrete gradient method (abbreviated to IA).

We set the starting point $\vec{x}^{0}=\vec{x}^{\delta}$, and the parameters to $\tau_{k}=1$ for all $k$, $\gamma=0.5$, and $\varphi_{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^{\delta}_{i}$, we use the scalar root solver scipy.optimize.brenth on Python. Otherwise, the updates are in closed form.

See Figure 7.5 for numerical results. By gradient norm, we mean $\operatorname{dist}(\partial F(\vec{x}^{k}),0)$. We observe that as in the convex case, the Bregman Itoh–Abe method converges significantly faster, due to the sparse structure of the problem.

In all of these cases, the sparsity structure when utilised properly leads to significantly faster convergence rates with the BSOR method.

8 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. 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 (lor14sfp, ). Furthermore, we will study corresponding inverse scale space schemes using other discrete gradients, such as the mean value discrete gradient.

Appendix A Counterexample of Theorem 4.5 for (3.7)

In the following example, we describe an example of a constrained optimisation problem for which the iterates of the 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\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{2}\to\mathbb{R}$ as described by Curry (cur44, ), for which the trajectory of a Euclidean gradient flow spirals along a gully, asymptotically trending towards the unit circle $S^{1}=\{[x,y]^{T}\in\mathbb{R}^{2}\;\mathrel{\mathop{\mathchar 58\relax}}\;x^{2}+y^{2}=1\}$.

Denote by $[\tilde{x}^{k},\tilde{y}^{k}]^{T}$ the iterates from the (standard) Itoh–Abe method for $W$ with (arbitrary) time steps $\tau\equiv 1$ and starting point $[\tilde{x}^{0},\tilde{y}^{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\subset\mathbb{R}^{2}$ with nonempty interior, such that $B\cap S^{1}\neq\emptyset$, and such that most of the iterates $([\tilde{x}^{k},\tilde{y}^{k}]^{T})_{k\in\mathbb{N}}$ are in $A$. That is, for all $k\in\mathbb{N}$, the set $\{i\leq k\;\mathrel{\mathop{\mathchar 58\relax}}\;[\tilde{x}^{i},\tilde{y}^{i}]^{T}\in A\}$ is strictly greater than the set $\{i\leq k\;\mathrel{\mathop{\mathchar 58\relax}}\;[\tilde{x}^{i},\tilde{y}^{i}]^{T}\notin 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\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{3}\to\mathbb{R}$ as folllows. Let $f_{A},f_{B}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{2}\to\mathbb{R}$ be smooth support functions of $A$ and $B$ respectively, with disjoint support sets. That is, $f_{A}(x,y)\in[0,1]$, $f_{B}(x,y)\in[0,1]$, and $f_{A}(x,y)\cdot f_{B}(x,y)=0$ for all $[x,y]^{T}\in\mathbb{R}^{2}$, and furthermore, $f_{A}(A)\equiv 1$ and $f_{B}(B)\equiv 1$. We then consider the optimisation problem

[TABLE]

We implement the unmofidied Bregman Itoh–Abe method (3.7) with $\tau\equiv 1$, $\vec{x}^{0}=[\tilde{x}^{0},\tilde{y}^{0},0]^{T}$, and $J(\vec{x})=\|\vec{x}\|^{2}/2+\chi_{[z\geq 0]}(z)$, where we use the shorthand notation $\vec{x}=[x,y,z]^{T}$.

If $[x,y]^{T}\in B\cap S^{1}$, then $\partial_{z}V(x,y,0)=-1$, meaning that $[x,y,0]^{T}$ is a nonstationary point of $V$ restricted to $\{z\geq 0\}$. Since the limit set of $([\tilde{x}^{k},\tilde{y}^{k}]])_{k\in\mathbb{N}}$ is $S^{1}$, it follows that the iterates $([\tilde{x}^{k},\tilde{y}^{k},0])_{k\in\mathbb{N}}$ admit nonstationary accumulation points for $V$. It therefore remains to show that for all $k\in\mathbb{N}$, $\vec{x}^{k}=[\tilde{x}^{k},\tilde{y}^{k},0]$ are admissible updates to the Bregman Itoh–Abe method (3.7) for (A.1).

To verify that $\vec{x}^{k}=[\tilde{x}^{k},\tilde{y}^{k},0]$ for all $k$, we can argue by induction. If $\vec{x}^{k}=[\tilde{x}^{k},\tilde{y^{k}},0]$, then $V(\tilde{x}^{k},\tilde{y}^{k},0)=W(\tilde{x}^{k},\tilde{y}^{k})$ so $\tilde{x}^{k+1}$ and $\tilde{y}^{k+1}$ are admissible updates for the $x$- and $y$-coordinates.

It remains to verify the same for the $z$-coordinate as well. If $[\tilde{x}^{k+1},\tilde{y}^{k+1}]^{T}\in A$, then $\partial_{z}V(\tilde{x}^{k+1},\tilde{y}^{k+1},0)=1$, so we have the update $z^{k+1}=0$ and $p^{k+1}_{3}=p^{k}_{3}-1\in\chi_{[z\geq 0]}(z^{k+1})$. On the other hand, if $[\tilde{x}^{k+1},\tilde{y}^{k+1}]^{T}\notin A$, then $\partial_{z}V(\tilde{x}^{k+1},\tilde{y}^{k+1},0)\in[-1,1]$, so if in addition $p^{k}_{3}<-1$, then we have the update $z^{k+1}=0$, $p^{k+1}_{3}=p^{k}_{3}-\partial_{z}V(\tilde{x}^{k+1},\tilde{y}^{k+1},0)$. Denoting by $M$ and $N$ the cardinalities of $\{i\leq k\;\mathrel{\mathop{\mathchar 58\relax}}\;[\tilde{x}^{i},\tilde{y}^{i}]^{T}\in A\}$ and $\{i\leq k\;\mathrel{\mathop{\mathchar 58\relax}}\;[\tilde{x}^{i},\tilde{y}^{i}]^{T}\notin A\}$ respectively, we know by construction of $A$ that $M-N\geq 1$. Furthermore, we observe that $p^{k}_{3}\leq N-M\leq-1$, which concludes the argument.

Appendix B $\ell^{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 $\tilde{x}^{k+1}_{i}$ the standard SOR update (5.2) for the $i$th coordinate. Then the $\ell^{1}$-regularised sparse SOR method can be expressed as follows.

If $x^{k}_{i}=0$ and $|\tilde{x}^{k+1}_{i}-\gamma r^{k}_{i}|\leq\gamma+\lambda\tau/a^{i}_{i}$, then

[TABLE] 2. 2.

Else if

[TABLE]

then

[TABLE] 3. 3.

Else if $x^{k}_{i}\neq 0$ and

[TABLE]

then set

[TABLE] 4. 4.

Else if $x_{i}\neq 0$ and

[TABLE]

then set

[TABLE]