∎\UseRawInputEncoding ¹²³⁴

Stochastic Bregman Subgradient Methods for Nonsmooth Nonconvex Optimization Problems

Kuangyu Ding

Kim-Chuan Toh (Received: date / Accepted: date)

Abstract

This paper focuses on the problem of minimizing a locally Lipschitz continuous function. Motivated by the effectiveness of Bregman gradient methods in training nonsmooth deep neural networks and the recent progress in stochastic subgradient methods for nonsmooth nonconvex optimization problems bolte2021conservative ; bolte2022subgradient ; xiao2023adam , we investigate the long-term behavior of stochastic Bregman subgradient methods in such context, especially when the objective function lacks Clarke regularity. We begin by exploring a general framework for Bregman-type methods, establishing their convergence by a differential inclusion approach. For practical applications, we develop a stochastic Bregman subgradient method that allows the subproblems to be solved inexactly. Furthermore, we demonstrate how a single timescale momentum can be integrated into the Bregman subgradient method with slight modifications to the momentum update. Additionally, we introduce a Bregman proximal subgradient method for solving composite optimization problems possibly with constraints, whose convergence can be guaranteed based on the general framework. Numerical experiments on training nonsmooth neural networks are conducted to validate the effectiveness of our proposed methods.

Keywords:

Nonsmooth nonconvex optimization Clarke regularity Conservative field Bregman subgradient methods Deep learning.

MSC:

90C05 90C06 90C25

Communicated by Sorin-Mihai Grad.

1 Introduction

In this paper, we focus on exploring the application of the Bregman distance based stochastic subgradient methods for solving nonsmooth nonconvex optimization problems. We consider the following unconstrained optimization problem,

$$\min_{x\in\mathbb{R}^{n}}f(x),$$

(1)

where f:ℝⁿ →ℝ is a locally Lipschitz continuous function (possibly nonconvex and nonsmooth). This class of problems encompasses many important applications, particularly in training nonsmooth neural networks where nonsmooth activate functions, such as the rectified linear unit (ReLU), are employed.

First order methods are commonly used to solve (1). In many scenarios, only noisy gradients or subgradients are available. During the last several decades, a great number of stochastic first order methods have been proposed. Among these, stochastic gradient descent (SGD) could be the most fundamental method. Based on SGD, a great number of its variants have also been developed to obtain some benefits on speed-up, stability, or memory efficiency. For example, the heavy-ball SGD polyak1964some , signSGD bernstein2018signsgd , and normalized SGD you2017large ; you2019large ; cutkosky2020momentum . These methods have demonstrated great efficiency and competitive generalization performance in various tasks in deep learning.

Despite the progress in developing SGD-type methods, the conventional convergence analysis largely pertains to scenarios where the objective function f exhibits certain regularity properties, namely, differentiability or weakly convexity. However, the prevalence of nonsmooth activation functions, such as ReLU or Leaky ReLU, in neural network architectures results in loss functions that often lack Clarke regularity (e.g. differentiability, weak convexity). Consequently, the conventional convergence analysis of SGD-type methods is not applicable in the context of training nonsmooth neural networks. Towards this issue, Bolte and Pauwels bolte2021conservative introduce the concept of conservative fields, which generalizes the Clarke subdifferential and admits a chain rule and sum rule for path-differentiable functions that may not be Clarke regular. Leveraging the concept of conservative fields and path-differentiability, subsequent studies bolte2021conservative ; castera2021inertial ; bolte2022subgradient ; le2023nonsmooth ; xiao2023adam ; xiao2023convergence ; ding2023adam ; xiao2024developing ; gurbuzbalaban2022stochastic ; ruszczynski2021stochastic utilize a differential inclusion approach benaim2005stochastic ; borkar2009stochastic to establish convergence for SGD-type methods in training nonsmooth deep neural networks. In particular, bolte2021conservative ; bolte2022subgradient study the convergence properties of SGD and proximal SGD for minimizing nonsmooth, path-differentiable functions. Moreover, castera2021inertial proposes the inertial Newton algorithm (INNA), which is regarded as a variant of momentum-accelerated SGD methods. Additionally, le2023nonsmooth establishes the convergence of SGD with heavy-ball momentum, and xiao2023convergence establishes the global stability of SGD with heavy-ball momentum. In addition, gurbuzbalaban2022stochastic ; ruszczynski2021stochastic design stochastic subgradient methods for solving multi-level nested optimization problems, and xiao2023adam ; ding2023adam establish the convergence properties of adaptive subgradient methods widely used in deep learning. Moreover, xiao2024developing develops subgradient methods to solve nonsmooth nonconvex constrained optimization problems. In summary, the crucial step in establishing convergence is to construct the differential inclusion corresponding to the iterative optimization algorithm and then to find a suitable Lyapunov function associated with that differential inclusion.

Beyond SGD-type methods, stochastic Bregman gradient methods, also known as mirror descent methods, have gained increasing interest recently. Initially introduced by Nemirovski and Yudin nemirovskij1983problem for solving constrained convex problems, the fundamental update scheme of Bregman gradient method is given as follows,

$$x_{k+1}=\mathop{\mathrm{argmin}}_{x\in\mathbb{R}^{n}}\left\lbrace \left\langle g_{k},x-x_{k}\right\rangle+\frac{1}{\eta_{k}}\mathcal{D}_{\phi}(x,x_{k})\right\rbrace ,$$

(2)

where D_ϕ represents the Bregman distance, a generalized measure of distance provided by a kernel function ϕ, with its definition to be presented in Section 2. A distinctive feature of Bregman gradient methods is the utilization of a broader selection of kernel functions, extending beyond the classical kernel function of 1/2‖ ·‖ ² used in SGD. Therefore, SGD is a special instance of Bregman gradient methods. The application of Bregman gradient methods has expanded to encompass both convex and nonconvex problems, as well as deterministic and stochastic contexts, as demonstrated by subsequent works beck2003mirror ; bauschke2017descent ; lu2018relatively ; bolte2018first ; zhang2018convergence ; dragomir2021optimal ; hanzely2021fastest ; yang2022bregman ; ding2023nonconvex . Bregman gradient methods have found many applications in optimization problems associated with probability constraints, such as optimal transport problem benamou2015iterative ; chu2023efficient and reinforcement learning lan2023policy ; zhan2023policy . Recent works have increasingly highlighted the potential of Bregman gradient methods in neural network training. Notably, ding2023nonconvex has shown that Bregman gradient methods can achieve better generalization performance and enhanced robustness in stepsize tuning compared to SGD for certain deep learning tasks. Additionally, many recent works, such as gunasekar2018characterizing ; wu2020continuous ; li2021implicit ; azizan2021stochastic ; sun2023unified , have demonstrated that Bregman gradient methods exhibit implicit regularization, leading to improved generalization performance with the selection of an appropriate kernel function. Additionally, amid2020reparameterizing ; ghai2022non ; li2022implicit provide an interpretation of the classical gradient descent method on reparameterized models from the perspective of Bregman gradient methods. This perspective has led to some applications in deep learning, including neural network quantization, as illustrated by ajanthan2021mirror .

However, existing convergence analysis for Bregman gradient methods has been limited to objective functions that are either differentiable or weakly convex. Given the increasing interest in applying Bregman gradient methods to train nonsmooth neural networks, and the current limited understanding of the convergence properties of Bregman subgradient methods for solving (1), this paper aims to provide a theoretical convergence guarantee for Bregman subgradient methods applied to nonsmooth nonconvex problems, particularly for the training of nonsmooth neural networks. Moreover, we propose practical Bregman subgradient methods specifically for training nonsmooth neural networks while ensuring their convergence.

To summarize, the contributions of this paper are as follows:

• General Bregman differential inclusion: We first investigate a general Bregman differential inclusion, whose discrete update scheme is consistent with that of Bregman-type methods and allows for biased evaluations of the abstract set-valued mapping. We establish the convergence for the discrete update scheme associated with this Bregman differential inclusion. Specifically, we demonstrate that any cluster point of the sequence generated by the discrete update scheme lies in the stable set of the Bregman differential inclusion, and the Lyapunov function values converge. A key aspect of our approach is the utilization of linear interpolation of the dual sequence induced by the kernel function.
• Applications of the general Bregman differential inclusion: By exploiting the flexibility of choosing the kernel function and the set valued mapping in the general framework, we introduce three types of stochastic Bregman subgradient methods for different scenarios. First, we consider the vanilla stochastic Bregman subgradient method (SBG) for unconstrained optimization problem, and establish its convergence. We further show that under certain regularity conditions of the kernel function, a preconditioned subgradient method can be regarded as an inexact Bregman subgradient method, thus fitting it within our proposed method. Moreover, we propose a single timescale momentum based stochastic Bregman subgradient method (MSBG), and establish its convergence by the sophisticated choice of the kernel function and set-valued mapping. Lastly, we consider the stochastic Bregman proximal subgradient method (SBPG) for constrained composite optimization problem and establish its convergence.
• Numerical experiments: To evaluate the performance of our proposed stochastic Bregman subgradient methods, we employ a block-wise polynomial kernel function in our proposed Bregman subgradient methods. We conduct numerical experiments to compare SGD, the momentum based stochastic Bregman subgradient method (MSBG), and its inexact version (iMSBG), for training nonsmooth neural networks. The results illustrate that our Bregman subgradient methods can achieve generalization performance comparable to that of SGD as well as the enhanced robustness of stepsize tuning.

The remaining sections of the paper are structured as follows. In Section 2, we provide preliminary materials on notations, nonsmooth analysis, and Bregman proximal mapping. Section 3 is dedicated to establishing the convergence of the general Bregman-type method based on the associated differential inclusion. The specific applications of the general Bregman-type method are presented in Section 4, where we first consider the vanilla stochastic Bregman subgradient method and then propose a single timescale momentum based Bregman subgradient method. The latter part of Section 4 explores the stochastic Bregman proximal subgradient method. Section 5 conducts numerical experiments to demonstrate the performance of our proposed Bregman subgradient methods. In the last section, we give a conclusion on this paper.

2 Preliminary

2.1 Notations

Given a proper and lower semicontinuous function f:ℝⁿ →ℝ̅:=(-∞,∞], we denote its domain as dom f={x:f(x)<∞}. The Fenchel conjugate function of f is defined as f∗(y):=sup{⟨x,y⟩-f(x):x∈ℝⁿ}. For a set S⊂ℝⁿ, we use cl S to denote its closure, int S to denote the set of its interior points, and conv(S) to denote its convex hull. A function f:S→ℝ is said to be of class Cᵏ(S) if it is k times differentiable and the k-th derivative is continuous on S. When there is no ambiguity regarding the domain, we simply use the notation Cᵏ. We let C(A,B) be the set of continuous mappings from set A to set B. We use ‖ ·‖ to denote the Euclidean norm for vectors and the Frobenius matrix norm for matrices. The d-dimensional unit ball is denoted by Bᵈ. The distance between a point w and a set A is denoted by dist (w,A):= ‖ w-u‖ : u∈ A. We use the convention dist(w,∅):=∞. The Minkowski sum of two sets A and B is denoted as A+B:={u+v: u∈ A, v∈ B}. We use α A:={α u: u∈ A} to denote the scaled set of A by a given scalar α. For a positive sequence {ηₖ}, we define λ_η(0):=0, λ_η(k):=Σᵢ₌₀ᵏ⁻¹ ηᵢ for k≥1, and Λ_η(t):=sup{k:λ_η(k)≤t}. In other word, Λ_η(t)=k if and only if λ_η(k)≤t<λ_η(k+1). Given a convex function g, we say x is a ν-optimal solution of the problem min_x∈ R ^n g(x) if dist(0,∂ g(x))≤ν.

Let (Ω,F,ℙ) be a probability space. Consider a stochastic process {ξₖ}ₖ≥0 and a filtration {Fₖ}ₖ≥0, where Fₖ is defined by the σ-algebra Fₖ:=σ(ξ₀,…,ξₖ) on Ω, the conditional expectation is denoted as 𝔼[·|Fₖ].

2.2 Nonsmooth analysis

In this subsection, we present some concepts from nonsmooth analysis, mainly based on clarke1990optimization ; RockWets98 .

Definition 1

Let X⊂ℝⁿ. A set-valued mapping S: X R ^m is said to be closed if its graph, defined by

$$\mathrm{graph}(S):=\lbrace (x,y)\in\mathbb{R}^{n}\times\mathbb{R}^{m}:y\in S(x)\rbrace$$

is a closed set in ℝⁿ ×ℝᵐ. Give a nonnegative δ, the δ-perturbed set of S is defined by

$$S^{\delta}(x):=\cup_{\lbrace y\in\mathbb{R}^{n}:\left\Vert y-x\right\Vert \leq\delta\rbrace }\left(S(y)+\delta\mathbb{B}^{m}\right).$$

S: X R ^m is called outer semicontinuous at x̅∈X if for any sequence x_i X →x and vᵢ ∈ S(xᵢ) converging to some v̅∈ℝᵐ, we have v̅∈ S(x̅). S is said to be outer semicontinuous if it is outer semicontinuous everywhere over X.

Definition 2

Let X⊂ℝⁿ be a closed set, the regular normal cone at x̅∈X is defined as N̂_X(x̅):={v∈ℝⁿ:⟨v,x-x̅⟩≤o(‖ x-x̅‖),for x∈X}. The limiting normal cone is defined by N^L_ X (x):= v∈ R ^n:x_k X →x,v_k∈N_ X (x_k),v_k→ v. The normal cone is defined by N_X(x̅):=cl conv(N^L_X(x̅)).

According to (RockWets98, , Proposition 6.6), the limiting normal cone N^L_X is outer semicontinuous, while this is not necessarily true for N_X.

Definition 3

Let f:ℝⁿ →ℝ be a locally Lipschitz continuous function. The generalized directional derivative of f at x∈ℝⁿ along the direction d∈ℝⁿ is defined by

$$f^{\circ}(x;d):=\lim\sup_{y\rightarrow x,t\downarrow 0}\frac{f(y+td)-f(y)}{t}.$$

The Clarke subdifferential of f at x is defined by

$$\partial f(x):=\left\lbrace v\in\mathbb{R}^{n}:\langle v,\,d\rangle\leq f^{\circ}(x;d),\;for\;all\;d\in\mathbb{R}^{n}\right\rbrace .$$

f is Clarke regular, if for any d∈ℝⁿ, its one-side directional derivative, defined by

$$f^{\ast }(x;d):=\lim_{t\downarrow 0}\frac{f(x+td)-f(x)}{t},$$

exists and f^(x;d)=f^ (x;d).

Clarke regularity excludes functions whose graph has upwards corners, such as -|x|. Notably, some basic calculus rules, including the sum and chain rules, may fail for the Clarke subdifferential without Clarke regularity. This limitation motivates the introduction of conservative fields.

2.2.1 Conservative field, path differentiability

In this part, we briefly introduce some relevant materials on conservative field, which are mainly based on bolte2021conservative .

Definition 4 (Conservative field and path-differentiability)

Let f:ℝⁿ →ℝ be a locally Lipschitz function. A set-valued mapping D_f: R ^n R ^n is a conservative field for f if it is nonempty closed valued, and has closed graph. For any absolutely continuous curve γ:[0,1]→ℝⁿ, f admits a chain rule with respect to D_f, i.e.

$$\frac{d}{dt}(f\circ\gamma)(t)=\langle v,\,\dot{\gamma}(t)\rangle,\text{for all}v\in D_{f}(\gamma(t))\text{and almost all}t\in[0,1].$$

(3)

If a locally Lipschitz function f admits a conservative field D_f, then we say that f is path-differentiable, and f is the potential function of D_f.

More generally, we can define conservative mappings for vector-valued functions, which serve as a generalization of the Jacobian.

Definition 5

Let F:ℝⁿ →ℝᵐ be a locally Lipschitz continuous function. J_F: R ^n R ^m× n is called a conservative mapping for F, if for any absolutely continuous curve γ:[0,1]→ℝⁿ, the following chain rule holds for almost all t∈[0,1]:

$$\frac{d}{dt}F(\gamma(t))=V\dot{\gamma}(t),\;for\;any\;V\in J_{F}(\gamma(t)).$$

Definition 6

Given a closed set X⊂ℝⁿ, we say that X admits a chain rule, if for any absolutely continuous curve γ:[0,1]→X, it holds that

$$\left\langle N_{\mathcal{X}}(\gamma(t)),\dot{\gamma}(t)\right\rangle=\lbrace 0\rbrace ,\text{for almost all}t\in[0,1].$$

As demonstrated in davis2020stochastic , when X is Whitney stratifiable (e.g. X is definable), then X admits the chain rule. We also use (d)/(dt)(fγ)(t)=⟨ D_f(γ(t)), γ(t)⟩ for almost all t∈[0,1] to represent that the chain rule (3) is valid for almost all t∈[0,1]. As stated in (bolte2021conservative, , Remark 3), if D_f is a conservative field for f, then conv(D_f) is also a conservative field for f, and D_f is locally bounded. Consequently, it follows from (RockWets98, , Theorem 5.7, Proposition 5.12) that for any x∈ℝⁿ and ε>0, there exist a neighborhood V of x, such that ∪_y∈ V D_f(y)⊂ D_f(x)+εBⁿ. It is important to note that the conservative field for a function f is not unique. Lemma 1 clarifies the relationship between a conservative field and the Clarke subdifferential, showing that the Clarke subdifferential is the smallest conservative field among all convex-valued conservative fields. We define x to be a D_f-stationary point if 0∈ D_f(x) to distinguish it from the conventional stationary point defined via the Clarke subdifferential. According to Lemma 1, a D_f-stationary point represents a weaker form of stationarity compared to a Clarke-stationary point. This weaker notion is inevitable if we want to capture the ”subgradients” produced by automatic differentiation (AD) in deep learning, for which the Clarke subdifferential may not capture in the absence of Clarke regularity. The same stationarity measure is also used in bolte2022subgradient ; le2023nonsmooth ; xiao2023adam .

Additionally, bolte2021conservative highlights the importance of the conservative field generated by the AD algorithm, widely utilized in modern deep learning frameworks like TensorFlow, PyTorch, and JAX. In the context when the conservative field corresponds precisely to the AD-generated subdifferential, we refer to such a D_f-stationary point as a stationary point in the sense of AD, which is the typical case in deep learning. To illustrate the conservative field generated by the AD algorithm, consider the loss function of a two-layer ReLU neural network:

$$f(W):=\frac{1}{2}\left\Vert W_{2}\sigma(W_{1}x)\right\Vert ^{2},$$

where σ(a):=max{a,0} denotes the (coordinate-wise) ReLU activation function, x is fixed, and W=(W₁,W₂). In Pytorch, the subgradient computed by AD corresponds to σ′(a)={(1,, a>0,; 0,, a≤0.). Consequently, if D_f represents the AD-generated conservative field, we have

$$D_{f}(W_{1},W_{2})=\left(\left(\sigma^{\prime}(W_{1}x)\circ(W_{2}^{T}W_{2}\sigma(W_{1}x))\right)x^{T},\;W_{2}\sigma(W_{1}x)\left(\sigma(W_{1}x)\right)^{T}\right).$$

Lemma 1

((bolte2021conservative, , Theorem 1, Corollary 1)) Consider a locally Lipschitz continuous function f:ℝⁿ →ℝ, with ∂ f as its Clarke subdifferential and D_f as its convex-valued conservative field. Then,

1. D_f(x)={∇ f(x)} for almost all x∈ℝⁿ.
2. ∂ f is a conservative field of f. Moreover, for any x∈ℝⁿ, it holds that ∂ f(x)⊂ D_f(x).

The next lemma provides the motivation for introducing the concept of conservative field, as it highlights that conservative fields preserve some basic calculus rules that do not necessarily hold for Clarke subdifferential without Clarke regularity.

Lemma 2

((bolte2021conservative, , Lemma 5, Corollary 4))

1. Let F₁:ℝⁿ →ℝᵐ, F₂:ℝᵐ →ℝˡ be locally Lipschitz continuous vector-valued functions, and J_1: R ^n R ^m× n and J_2: R ^m R ^l× m be the conservative mappings of F₁ and F₂ respectively. Then the mapping J_2 J_1: R ^n R ^l× n is a conservative mapping for F_2 F_1.
2. Let D_fᵢ be a conservative field for fᵢ, i=1,...,N. Then, D_f=Σᵢ₌₁^N D_fᵢ is a conservative field of f=Σᵢ₌₁^N fᵢ.

The Clarke subdifferential may fail to satisfy key calculus rules, whereas the conservative field retains these properties, as established in Lemma 2. This ensures the well-definedness of gradient methods based on automatic differentiation (AD) in deep learning，since ”subdifferential” generated by AD algorithms is a special case of conservative field. A notable subclass of path-differentiable functions is the set of definable functions in o-minimal structures coste2000introduction ; van1996geometric , which includes widely used activation and loss functions, such as sigmoid, ReLU, softmax, l₁ loss, cross-entropy, and hinge loss. Definable functions are closed under composition, summation, and set-valued integration, ensuring that neural networks built from such blocks remain definable davis2020stochastic ; bolte2021conservative . Their geometric structure, characterized by Whitney stratification, enables effective analysis of nonsmooth optimization algorithms bolte2007clarke ; davis2020stochastic ; bolte2021conservative .

2.3 Bregman proximal mapping

In this subsection, we review some concepts related to Bregman proximal mappings. Given that this paper mainly focuses on unconstrained problems, we restrict our discussion on the entire space ℝⁿ. For more general concepts about Bregman gradient mapping, readers can refer to works such as bauschke2017descent ; bolte2018first ; ding2023nonconvex .

Definition 7 (Kernel function and Bregman distance over ℝⁿ). A function ϕ:ℝⁿ →ℝ is called a kernel function over ℝⁿ, if ϕ is convex and ϕ∈C¹(ℝⁿ). The Bregman distance bregman1967relaxation generated by ϕ is denoted as D_ϕ(x,y):ℝⁿ ×ℝⁿ →[0,+∞), where

$$\mathcal{D}_{\phi}(x,y)=\phi(x)-\phi(y)-\langle\nabla\phi(y),\,x-y\rangle.$$

Definition 8 (Legendre kernel over ℝⁿ). Let ϕ be a kernel function over ℝⁿ, such that limₖ→∞ ‖ ∇ϕ(xₖ)‖=∞ whenever {xₖ}ₖ∈ℕ satisfies limₖ→∞ ‖ xₖ ‖=∞. The function ϕ is called a Legendre function over ℝⁿ if it is also strictly convex on ℝⁿ.

Definition 9

Given a locally Lipschitz continuous function R and a Legendre kernel function ϕ∈C¹(ℝⁿ), we denote the Bregman proximal mapping by Prox^ϕ_R:=(∇ϕ+∂ R)⁻¹ ∇ϕ, which is a set-valued mapping defined as follows,

$$\mathrm{Prox}^{\phi}_{R}(x):=\mathop{\mathrm{argmin}}_{u\in\mathbb{R}^{n}}\;\lbrace R(u)+\mathcal{D}_{\phi}(u,x)\rbrace .$$

(4)

Under mild conditions, Prox^ϕ_R(x) is a nonempty compact set for any x∈ℝⁿ, which will ensure the well-posedness of our methods. We have the following lemma, which directly follows from Weierstrass’s theorem, we omit its proof for simplicity.

Lemma 3

Let R be a continuous function and ϕ be a Legendre function over ℝⁿ. Suppose α R+ϕ is supercoercive, i.e. lim_‖ x‖ →∞(α R(x)+ϕ(x))/(‖ x‖)=∞, for any α>0. Then, for any x∈ℝⁿ, the set Prox^ϕ_R(x) is a nonempty compact subset of ℝⁿ.

In the context of Bregman proximal gradient method, we focus on the Bregman forward-backward splitting operator T ^ϕ_α,R: R ^n× R ^n R ^n defined by:

$$\mathrm{T}^{\phi}_{\alpha,R}(x,v):=\mathop{\mathrm{argmin}}_{z\in\mathbb{R}^{n}}\;\left\lbrace \langle v,\,z-x\rangle+R(z)+\frac{1}{\alpha}\mathcal{D}_{\phi}(z,x)\right\rbrace .$$

For any given v∈ℝⁿ and α>0, under the same assumption in Lemma 3, T^ϕ_α,R is also well defined.

We end this subsection by introducing the concept of functional convergence over any compact set, as discussed in literature such as benaim2005stochastic ; duchi2018stochastic ; davis2020stochastic . Given a sequence of mappings fₙ ∈C(ℝ₊,ℝⁿ), we say that fₙ converges to f in C(ℝ₊,ℝⁿ) if, for any T>0, it holds that

$$\lim_{n\rightarrow\infty}\sup_{t\in[0,T]}\Vert f_{n}(t)-f(t)\Vert =0.$$

3 Bregman differential inclusion

In this section, we investigate the Bregman-type differential inclusion and its discrete approximation, which coincides with the iterative Bregman subgradient methods. The analysis tools employed in this section draw inspiration from various works on stochastic approximation, such as benaim2005stochastic ; borkar2009stochastic ; duchi2018stochastic ; davis2020stochastic ; ruszczynski2020convergence . To establish the convergence of the discrete sequence {xₖ}, a key idea is to show that the linear interpolation of the sequence {xₖ} defined by

$${z(t)}:=x_{k}+\frac{t-\lambda_{\eta}(k)}{\eta_{k}}(x_{k+1}-x_{k}),\;t\in[\lambda_{\eta}(k),\lambda_{\eta}(k+1))$$

(5)

is a perturbed solution benaim2005stochastic to the associated differential inclusion, where {ηₖ} serves as the stepsize in the subgradient methods. However, due to the non-Euclidean nature of Bregman subgradient methods, this methodology requires modification. Given a kernel function ϕ and a general set-valued mapping H, we consider the following general differential inclusion:

$$\frac{d}{dt}\nabla\phi(x(t))\in-\mathcal{H}(x(t)),\text{for almost all}t\geq 0\text{.}$$

(6)

Any absolutely continuous solution to (6) is termed a trajectory of (6). The stable set of (6) is defined as

$$\mathcal{H}^{-1}(0):=\lbrace x\in\mathbb{R}^{n}:0\in\mathcal{H}(x)\rbrace .$$

(7)

For any C² convex function ϕ, the differential inclusion (6) can be interpreted as a gradient flow equipped with the Riemannian metric induced by ⟨·,·⟩_∇² ϕ(x), as demonstrated in works such as bolte2003barrier ; alvarez2004hessian . The corresponding discrete scheme of (6) is given by

$$\nabla\phi(x_{k+1})=\nabla\phi(x_{k})-\eta_{k}(d_{k}+\xi_{k}),$$

(8)

where dₖ is an evaluation of H(xₖ) with possible inexactness, and ξₖ is the stochastic noise. This formulation, referred to as the general Bregman-type method, is notable for the versatile choices of ϕ and H. For the Bregman counterpart of the interpolated process (5) in the Euclidean setting, we introduce the linear interpolation for the dual sequence {∇ϕ(xₖ)}:

$$x(t):=\nabla\phi^{\ast }\left(\nabla\phi(x_{k})+\frac{t-\lambda_{\eta}(k)}{\eta_{k}}(\nabla\phi(x_{k+1})-\nabla\phi(x_{k}))\right),\;t\in[\lambda_{\eta}(k),\lambda_{\eta}(k+1)).$$

(9)

If ϕ∗ ∈C¹(ℝⁿ), then (9) is well defined. Let xᵗ(·) denote the time-shifted curve of the interpolated process, i.e. xᵗ(·)=x(t+·). We make the following assumptions on (6) and (8) to ensure that the iterative sequence generated by (8) tracks a trajectory of (6) asymptotically.

Assumption 3.1

1. ϕ is a supercoecive Legendre kernel function over ℝⁿ, and ∇ϕ is differentiable almost everywhere.
2. The sequences {xₖ}, {∇ϕ(xₖ)} and {dₖ} are uniformly bounded.
3. The stepsize {ηₖ} satisfies Σₖ₌₀^∞ ηₖ=∞ and limₖ→∞ ηₖ=0.
4. For any T>0, the noise sequence {ξₖ} satisfies

$$\lim_{s\rightarrow\infty}\sup_{s\leq i\leq\Lambda_{\eta}(\lambda_{\eta}(s)+T)}\left\Vert \sum_{k=s}^{i}\eta_{k}\xi_{k}\right\Vert =0.$$

(10)

5. The set-valued mapping H has a closed graph. Additionally, for any unbounded increasing sequence {kⱼ} such that {x_kⱼ} converges to x̅, it holds that

$$\lim_{N\rightarrow\infty}\mathrm{dist}\left(\frac{1}{N}\sum_{j=1}^{N}d_{k_{j}},\mathcal{H}(\bar{x})\right)=0.$$

(11)

Here are some remarks on Assumption 3.1.

Remark 1

1. If the Legendre function ϕ is supercoecive, i.e. lim_‖ u‖ →∞(ϕ(u))/(‖ u‖)=∞, then by (rockafellar1997convex, , Theorem 26.5, Corollary 13.3.1), ϕ∗ ∈C¹(ℝⁿ) is strictly convex, and (∇ϕ)⁻¹=∇ϕ∗.

2. Uniform boundedness of either {xₖ} or {∇ϕ(xₖ)} may suffice under mild conditions. For example, if ϕ is locally strongly convex, then the uniform boundedness of {∇ϕ(xₖ)} implies the uniform boundedness of {xₖ}. Conversely, if ∇ϕ is locally Lipschitz continuous, then the uniform boundedness of {xₖ} leads to the uniform boundedness of {∇ϕ(xₖ)}.

3. For a martingale difference noise sequence {ξₖ}, i.e. 𝔼[ξₖ₊₁|Fₖ]=0 holds almost surely, as shown in benaim2005stochastic , (10) can be ensured almost surely under one of the following conditions:

   1. Uniform boundedness of {ξₖ} with ηₖ=o((1)/(log k)) as utilized in castera2021inertial ; xiao2023adam ; xiao2023convergence .
   2. Variance-bounded {ξₖ}, i.e. 𝔼[‖ ξₖ₊₁ ‖ ²|Fₖ]≤σ² for some σ>0, and square-summable {ηₖ}, i.e. Σₖ₌₀^∞ ηₖ²<∞, as employed in davis2020stochastic .

4. Assumption 3.1(5) describes how dₖ approximates H(xₖ). This assumption is quite mild. A sufficient condition for it hold is that dₖ ∈Hᵟₖ(xₖ) with limₖ→∞ δₖ=0.

The following theorem demonstrates that the discrete sequence generated by the general Bregman-type method (8) asymptotically tracks a trajectory of the differential inclusion (6). Our key technical novelty, compared to previous works on stochastic approximation, especially benaim2005stochastic , is in establishing convergence under weaker assumptions on the kernel regularity. In benaim2005stochastic , the analysis focuses on differential inclusions of the form ẋ∈-G(x) and discrete updates are explicitly defined by xₖ₊₁=xₖ-ηₖ gₖ with gₖ ∈G(xₖ). Directly applying their results to our Bregman-type updates would require verifying convergence to a set-valued mapping H̃, specifically establishing:

$$\mathrm{dist}\left(\frac{x_{k}-\nabla\phi^{\ast }(\nabla\phi(x_{k})-\eta_{k}d_{k})}{\eta_{k}},\widetilde{\mathcal{H}}(x_{k})\right)\to 0.$$

When the kernel ϕ is C²(ℝⁿ), this convergence can indeed be guaranteed, as we can explicitly choose H̃(x)=(∇² ϕ(x))⁻¹ H(x). However, our analysis relaxes this requirement significantly, assuming only that ∇ϕ is differentiable almost everywhere. Under this weaker condition, such a direct equivalence no longer holds. To overcome this issue, we adopt a carefully designed interpolation scheme for the dual sequence, as defined in (9). By employing this novel linear interpolation approach on the sequence {∇ϕ(xₖ)}, we successfully establish convergence without the need for twice continuous differentiability of ϕ. The complete details of our approach and proofs are provided in Appendix A.

Theorem 3.1

Suppose Assumption 3.1 holds. Then, for any sequence {τₖ}ₖ₌₁^∞ ⊂ℝ₊, the set of shifted sequence {∇ϕ(x^τₖ(·))}ₖ₌₁^∞ is relatively compact in C(ℝ₊,ℝⁿ). If limₖ→∞ τₖ=∞, then any cluster point y̅(·) of {∇ϕ(x^τₖ(·))}ₖ₌₁^∞ belongs to C(ℝ+,ℝⁿ). Define x̅(·):=∇ϕ∗(y̅(·)). Then there exists a measurable d̅(t)∈H(x̅(t)) satisfying

$$\nabla\phi(\bar{x}(t))=\nabla\phi(\bar{x}(0))-\int_{0}^{t}\bar{d}(\tau)d\tau\quad for\;all\;t\geq 0.$$

(12)

Equivalently, we have

$$\frac{d}{dt}\nabla\phi(\bar{x}(t))\in-\mathcal{H}(\bar{x}(t)),\;for\;almost\;all\;t\geq 0.$$

The following assumption ensures that the trajectory subsequentially converges to the stable set of (6), and the Lyapunov function values converge.

Assumption 3.2

There exists a continuous function Ψ:ℝⁿ →ℝ, such that the following conditions hold:

1. (Weak Morse-Sard). The set {Ψ(x):for x such that 0∈H(x)} has empty interior in ℝ.
2. (Lyapunov function). Ψ is lower bounded, i.e. _x∈ R ^n(x)> -∞. For any trajectory z(t) of the differential inclusion (6) with z(0)∉H⁻¹(0), there exists T>0 such that

Ψ(z(T))<sup_[t∈[0,T]] Ψ(z(t))≤Ψ(z(0)).

We say that a continuous function Ψ, satisfying Assumption 3.2.2, is a Lyapunov function for the differential inclusion (6) with a stable set H⁻¹(0).

Now, we are ready to present the main theorem in this section.

Theorem 3.2

Suppose Assumptions 3.1 and 3.2 hold. Then any cluster point of {xₖ} lies in H⁻¹(0) and the function values {Ψ(xₖ)}ₖ≥1 converge.

Proof

Since {xₖ} is bounded, let x∗ be any cluster point of {xₖ}, and limᵢ→∞x_kᵢ=x∗. By Theorem 3.1, up to a subsequence, ∇ϕ(x^λ_η(kᵢ)(·))→∇ϕ(x̅(·)) for some x̅(·)∈C(ℝ₊,ℝⁿ). Note that x_kᵢ=x^λ_η(kᵢ)(0), so we have

$${} \left\Vert \bar{x}(0)-x^{\ast }\right\Vert \leq \limsup_{i\rightarrow\infty}\left(\left\Vert \bar{x}(0)-x^{\lambda_{\eta}(k_{i})}(0)\right\Vert +\left\Vert x^{\lambda_{\eta}(k_{i})}(0)-x_{k_{i}}\right\Vert +\left\Vert x_{k_{i}}-x^{\ast }\right\Vert \right)\qquad {} = \limsup_{i\rightarrow\infty}\left(\left\Vert \nabla\phi^{\ast }(\nabla\phi(\bar{x}(0)))-\nabla\phi^{\ast }(\nabla\phi(x^{\lambda_{\eta}(k_{i})}(0)))\right\Vert +\left\Vert x^{\lambda_{\eta}(k_{i})}(0)-x_{k_{i}}\right\Vert +\left\Vert x_{k_{i}}-x^{\ast }\right\Vert \right)\qquad {} = 0,$$

where the last equality comes from the continuity of ∇ϕ∗. Hence x̅(0)=x∗. Suppose x∗∉H⁻¹(0), then by Assumption 3.2, there exists T>0, such that

Ψ(x̅(T))<sup_[t∈[0,T]] Ψ(x̅(t))≤Ψ(x∗).

On the other hand, by Lemma B.10, Ψ(x(t)) converges as t→∞. Therefore, we obtain

$$\Psi(\bar{x}(T))={\lim_{i\rightarrow\infty}}\Psi(x^{\lambda_{\eta}(k_{i})}(T))={\lim_{i\rightarrow\infty}}\Psi(x({\lambda_{\eta}(k_{i})}+T))=\lim_{t\rightarrow\infty}\Psi(x(t))=\Psi(x^{\ast }),$$

which leads to a contradiction. Therefore, 0∈H(x∗). This completes the proof.

The results in Theorem 3.2 illustrate that the sequence {xₖ} finds the stable set of (6), and the function values of Lyapunov function converge. Due to the versatility in selecting the set-valued mapping H, we will later show that if H is chosen as the conservative field of f, the sequence generated by the Bregman-type method (8) converges to a stationary point of the optimization problem (1).

In the remainder of this section, we aim to establish the global stability of the framework (8) under the noise sequence {ξₖ} correspond to random reshuffling. To that end, we impose the following assumptions on the stepsize {ηₖ} and the noise {ξₖ}.

Assumption 3.3

There exists an integer N>0 such that

1. For any nonnegative integers i,j<N, it holds that ηₖN+i=ηₖN+j for any k∈ℕ₊.
2. The sequence {ξₖ} is uniformly bounded. Moreover, for any j∈ℕ₊, almost surely, it holds that Σ_k=jN⁽j+1)N-1 ξₖ₊₁=0.

Here we make some remarks on Assumptions 3.3. Assumption 3.3(1) indicates that the stepsizes {ηₖ} remain constant within each epoch, which is a standard setting in neural network training tasks. Moreover, Assumption 3.3(2) encodes the essence of random reshuffling in the stochastic subgradient computation of {dₖ} in (8). Section 4 provides a concrete example in the context of finite-sum minimization, illustrating why this assumption is reasonable.

Lemma 4

Suppose Assumption 3.3 holds for the sequence of noises {ξₖ} and stepsizes {ηₖ}. Moreover, the Lyapunov function Ψ associated with H is coercive. Then for any ε>0 and T>0, there exists η_ε>0 such that for any {ηₖ} satisfying _k→+∞_k≤_ε, it holds that

$$\limsup_{s\to+\infty}\sup_{s\leq i\leq{\Lambda_{\eta}(\lambda_{\eta}(s)+T)}}\left\Vert \sum_{k=s}^{i}\eta_{k}\xi_{k+1}\right\Vert \leq\varepsilon.$$

(13)

We now state a theorem establishing the global stability of (8) under non-diminishing stepsizes {ηₖ}.

Theorem 3.3

Suppose Assumptions 3.2 and 3.3 hold. Additionally, assume ϕ is strongly convex and Lipschitz smooth, and that dₖ ∈Hᵟₖ(xₖ) with limₖ→∞ δₖ=0. Then for any ε>0, there exists ηₘₐₓ(ε)>0 such that for any {ηₖ} satisfying _k→+∞_k≤_max(ε), the following holds:

$$\limsup_{k\to+\infty}\;\mathrm{dist}\left(x_{k},\lbrace x\in\mathbb{R}^{n}:0\in\mathcal{H}(x)\rbrace \right)\leq\varepsilon.$$

(14)

Consequently, the sequence {xₖ} is uniformly bounded.

Theorem 3.3 implies that as long as the stepsizes {ηₖ} are sufficiently small, the sequence {xₖ} remains stable. Notably, when random reshuffling noise is present, together with other essential mild conditions, the uniform boundedness of the iterates can be rigorously established, thereby justifying the boundedness assumption in Assumption 3.1.

4 Applications

Based on the framework of the general Bregman-type method as outlined in (6), in this section, we consider three specific types of stochastic Bregman subgradient methods by choosing different types of kernel function ϕ and set-valued mapping H. In the first two parts, we consider vanilla and single timescale momentum based stochastic Bregman subgradient methods for unconstrained optimization problems. Subsequently, we extend our methods to the stochastic Bregman proximal subgradient method for solving constrained composite optimization problems.

4.1 Stochastic Bregman subgradient method

In this subsection, we consider the following stochastic Bregman subgradient update scheme:

$$x_{k+1} \approx\arg\min_{x\in\mathbb{R}^{n}}\;\left\lbrace \left\langle g_{k},x-x_{k}\right\rangle+\frac{1}{\eta_{k}}\mathcal{D}_{\phi}(x,x_{k})\right\rbrace ,\qquad x_{k+1}\quad \nu_{k}$$

where gₖ=dₖ+ξₖ, dₖ ∈ Dᵟₖ_f(xₖ), and ξₖ is the stochastic noise. The associated differential inclusion is

$$\frac{d}{dt}\nabla\phi(x(t))\in-D_{f}(x(t)).$$

(15)

Recall the definition of ν-optimal solution in Section 2.1, xₖ₊₁ is a νₖ-optimal solution of the subproblem, if and only if ‖ gₖ+(1)/(ηₖ)(∇ϕ(xₖ₊₁)-∇ϕ(xₖ))‖≤νₖ. Given the allowance for inexact solutions in the SBG framework, we illustrate that a kernel Hessian preconditioned subgradient method fits within our framework, akin to the approach in the recently proposed ABPG in takahashi2023approximate . The concept of Hessian preconditioning has also been examined in continuous settings as seen in alvarez2004hessian ; bolte2003barrier . Assuming the absence of the nonsmooth term in ABPG, the existence and nonsingularity of ∇² ϕ everywhere, and that supₖ ‖ (∇² ϕ(xₖ))⁻¹ ‖≤c for some c>0, the ABPG updates scheme in takahashi2023approximate is given by

$$x_{k}^{+}=x_{k}-\eta_{k}(\nabla^{2}\phi(x_{k}))^{-1}g_{k}.$$

(16)

This yields

$$\nabla\phi(x_{k}^{+})-(\nabla\phi(x_{k})-\eta_{k}g_{k})\qquad {} = \nabla\phi(x_{k}-\eta_{k}(\nabla^{2}\phi(x_{k}))^{-1}g_{k})-(\nabla\phi(x_{k})-\eta_{k}g_{k})\qquad {} = \nabla\phi(x_{k})-\eta_{k}\nabla^{2}\phi(x_{k})(\nabla^{2}\phi(x_{k}))^{-1}g_{k}+o(\eta_{k})-(\nabla\phi(x_{k})-\eta_{k}g_{k})\qquad {} = o(\eta_{k}),$$

indicating that limₖ→∞ ‖ gₖ+(∇ϕ(x⁺ₖ)-∇ϕ(xₖ))/(ηₖ)‖=0, which implies that x⁺ₖ is an approximate solution to the SBG subproblem.

We make the following assumptions on (SBG).

Assumption 4.1

1. ϕ is a supercoecive Legendre kernel function over ℝⁿ, and ∇ϕ is differentiable almost everywhere. Moreover, for any absolutely continuous mapping z(·)∈C(ℝ₊,ℝⁿ), ∇² ϕ(z(s)) is positive definite for almost all s≥0.

2. The sequences {xₖ}, {∇ϕ(xₖ)} and {dₖ} are uniformly bounded almost surely.

3. {ξₖ} is a martingale difference noise, i.e. 𝔼[ξₖ₊₁|Fₖ]=0 holds almost surely. Σₖ₌₀^∞ ηₖ=∞, and the stepsize and noise satisfy one of the following two conditions:

   1. {ξₖ} is uniformly bounded and ηₖ=o((1)/(log k)).
   2. {ξₖ} has bounded variance, i.e. 𝔼[‖ ξₖ₊₁ ‖ ²|Fₖ]≤σ²<∞, and Σₖ₌₀^∞ ηₖ²<∞.

4. limₖ→∞ δₖ=0, and limₖ→∞ νₖ=0.

Now, we make some remarks regarding Assumption 4.1. Assumption 4.1(1) and (3) are directly adapted from Assumption 3.1(1) and (4), respectively. Assumption 4.1(4) characterizes how dₖ approximates D_f(xₖ) and quantifies the inexactness in solving the subproblem. Assumption 4.1(2) can be ensured by global stability under randomly reshuffling noise {ξₖ}, as demonstrated by Theorem 3.3. We elaborate on these points further below.

Consider the case where f(x)=(1)/(N)Σᵢ₌₁^N fᵢ(x) being coercive, each fᵢ is path-differentiable and admits a convex-valued local bounded conservative field D_fᵢ(x). At each iteration, we sample an index iₖ from [N]:={1,2,…,N} such that for each epoch j≥0, {iₖ:jN≤k<(j+1)N}=[N] holds. We then compute dₖ ∈ D_f_iₖ(xₖ). This setup corresponds to the standard random reshuffling framework. For each k≥0, let jₖ≥0 be such that k∈[jₖ N,(jₖ+1)N), it follows that dₖ ∈ D_f(xₖ)+ξₖ₊₁, where ξₖ₊₁ represents the evaluation noise. Defining δₖ=2NΣ_jₖ N≤l<(jₖ+1)N‖ xₗ₊₁-xₗ ‖, we have dₖ ∈ D_f_iₖ(xₖ)⊆ D_f_iₖᵟₖ/(2N)(x_jₖ N). Consequently,

$$\frac{1}{N}\sum_{j_{k}N\leq l\lt (j_{k}+1)N}d_{l}\in\frac{1}{N}\sum_{j_{k}N\leq l\lt (j_{k}+1)N}D_{f_{i_{l}}}^{\delta_{k}/(2N)}(x_{j_{k}N})\subseteq D_{f}^{\delta_{k}/2}(x_{j_{k}N})\subseteq D_{f}^{\delta_{k}}(x_{k}).$$

Furthermore, it holds that (1)/(N)ΣⱼN≤k<(j+1)N ξₖ₊₁=0. Therefore, we can conclude that Assumption 3.3(2) corresponds to the setting where the stochastic subgradients {dₖ} are generated by the random reshuffling sampling technique. Therefore, applying Theorem 3.3, we can obtain the uniform boundedness of {xₖ}. Since D_f is locally bounded and ϕ is locally Lipschitz (implying that ∇ϕ is locally bounded), the boundedness of the sequence {xₖ} guarantees that both {∇ϕ(xₖ)} and {dₖ} are uniformly bounded (see, e.g., (RockWets98, , Proposition 5.15)).

To ensure the convergence of (SBG), we make the following assumptions on f.

Assumption 4.2

1. f is lower bounded, i.e. _x∈ R ^nf(x)> -∞. Moreover, f is a potential function that admits D_f as its convex-valued conservative field.
2. The critical value set {f(x): 0∈ D_f(x)} has empty interior in ℝ.

Assumption 4.2 is a specific instance of Assumption 3.2. In particular, f is the Lyapunov function and Assumption 4.2(2) requires a weak Sard condition for f. Recall that the classical Sard theorem for differentiable functions in ℝⁿ states that the set of critical values has Lebesgue measure zero.

We have the following two propositions.

Lemma 5

Suppose Assumptions 4.1 and 4.2 hold. For any dᵉₖ such that ‖ dᵉₖ ‖≤νₖ, and any increasing sequence {kⱼ} such that {x_kⱼ} converges to x̅, it holds that

$$\lim_{N\rightarrow\infty}\mathrm{dist}\left(\frac{1}{N}\sum_{j=1}^{N}(d_{k_{j}}+d^{e}_{k_{j}}),D_{f}(\bar{x})\right)=0.$$

Proof

By the inexact condition in (SBG), it follows that there exists dᵉₖ, such that ‖ dᵉₖ ‖≤νₖ and

$$\nabla\phi(x_{k+1})=\nabla\phi(x_{k})-\eta_{k}\left(d_{k}+d^{e}_{k}+\xi_{k}\right).$$

Define d̃ₖ:=dₖ+dᵉₖ ∈ D_fᵟ̃ₖ(xₖ), where δ̃ₖ=δₖ+νₖ. Note that limₖ→∞δ̃ₖ=0. Since D_f has a closed graph, then for any {x_kⱼ} converging to x̅, it holds that limⱼ→∞ dist(d̃_kⱼ,D_f(x̅))=0. Note that D_f(x̅) is a convex set, by Jensen’s inequality, we have

$$\lim_{N\rightarrow\infty}\mathrm{dist}\left(\frac{1}{N}\sum_{j=1}^{N}\tilde{d}_{k_{j}},D_{f}(\bar{x})\right)=\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{j=1}^{N}\mathrm{dist}\left(\tilde{d}_{k_{j}},D_{f}(\bar{x})\right)=0.$$

This completes the proof.

Proposition 1

Suppose Assumptions 4.1 and 4.2 hold. Then f is a Lyapunov function for the differential inclusion (15) with the stable set {x∈ℝⁿ: 0∈ D_f(x)}.

Proof

Consider any trajectory z(·) for the differential inclusion (15) with 0∉D_f(z(0)). We have that for almost all s≥0,

$$\frac{d}{ds}f(z(s))=\left\langle D_{f}(z(s)),\dot{z}(s)\right\rangle\ni-\left\langle\nabla^{2}\phi(z(s))\dot{z}(s),\dot{z}(s)\right\rangle$$

Therefore, for any t≥0, it holds that

$$f(z(t))-f(z(0))=-\int_{0}^{t}\left\langle\nabla^{2}\phi(z(s))\dot{z}(s),\dot{z}(s)\right\rangle ds\leq-\int_{0}^{t}\lambda_{\min}\left(\nabla^{2}\phi(z(s))\right)\left\Vert \dot{z}(s)\right\Vert ^{2}ds\leq 0.$$

We now prove the required result by contradiction. Suppose for any t≥0, f(z(t))=f(z(0)), then, we have λ_min(^2ϕ(z(s)))‖ z(s)‖ ^2=0 for almost all s≥0. Since λ_min(^2ϕ(z(·)))> 0 almost everywhere in ℝ₊, then ż(s)=0 for almost all s≥0. Since z(·) is absolutely continuous, therefore, z(t)≡z(0) for any t≥0. Then, 0=(d)/(dt)∇ϕ(z(t))∈-D_f(z(t))=-D_f(z(0)). This is contradictory to the fact that z(0) is not a D_f-stationary point of f. Therefore, there exists T>0, such that f(z(T))<sup_[t∈[0,T]] f(z(t))≤f(z(0)). This completes the proof.

By Lemma 5, Proposition 1 and Theorem 3.2, we can directly derive the following convergence results for (SBG).

Theorem 4.1

Suppose Assumptions 4.1 and 4.2 hold. Then almost surely, any cluster point of {xₖ} generated by (SBG) is a D_f-stationary point and the function values {f(xₖ)} converge.

4.2 Momentum based stochastic Bregman subgradient method

In this section, we introduce a momentum-based stochastic Bregman subgradient method (MSBG). The momentum technique is widely adopted in practice, particularly for training neural networks using stochastic first-order methods such as SGDM, Adam, and AdamW, primarily due to its empirical effectiveness. Nevertheless, a complete theoretical understanding of momentum’s benefits, even within smooth optimization contexts, remains elusive. While certain studies have attempted to interpret momentum through generalization perspectives such as margin analysis and algorithmic stability (jelassitowards ; ramezani2024generalization ), such analysis is beyond the scope of pure optimization, and thus we omit further discussion of these issues here.

From an optimization perspective, several studies (defazio2020understanding ; liu2020improved ; guo2021novel ) highlight that momentum can accelerate convergence by canceling out some impact of stochastic gradient noise. Intuitively, momentum reduces variance in the gradient estimates, effectively reducing the noise that often dominates convergence behavior, particularly in high-noise settings typical of neural network training. Motivated by these observations, we also incorporate momentum into our proposed Bregman subgradient method (MSBG), leveraging its variance-reducing benefits to potentially enhance convergence performance in practice.

For a chosen kernel function φ:ℝⁿ →ℝ, the momentum based update scheme is given as follows:

$$\left\lbrace x_{k+1}\approx\mathop{\mathrm{argmin}}_{x\in\mathbb{R}^{n}}\left\lbrace \left\langle m_{k},x-x_{k}\right\rangle+\frac{1}{\eta_{k}}\mathcal{D}_{\varphi}(x,x_{k})\right\rbrace\quad \text{s.t.}{x_{k+1}\text{is a}\nu_{k}\text{-optimal solution of the subproblem.}}\quad m_{k+1}=m_{k}-\theta_{k}P(x_{k})(m_{k}-g_{k}), \right.$$

(MSBG)

where gₖ=dₖ+ξₖ, dₖ ∈ Dᵟₖ_f(xₖ), and P(xₖ)∈ℝⁿ× n denotes a preconditioning matrix. Similar to (16), the MSBG subproblem can also be solved in an inexact manner by adopting a preconditioned subgradient strategy as shown in (16). For the ease of presentation, we omit the discussion. Our MSBG method is a single timescale method in the sense that the stepsize ηₖ for the primal variable and the stepsize θₖ for the momentum decay at the same rate.

We make the following assumptions on (MSBG).

Assumption 4.3

1. φ∈C²(ℝⁿ) is a supercoecive Legendre kernel function over ℝⁿ, and ∇² φ(·) is positive definite everywhere. Moreover, P(·)=(∇² φ(·))⁻¹:ℝⁿ →ℝⁿ× n.

2. The sequences {xₖ}, {∇ϕ(xₖ)}, {dₖ}, and {mₖ} are uniformly bounded almost surely.

3. {ξₖ} is a martingale difference noise, Σₖ₌₀^∞ ηₖ=∞, and the stepsize and noise satisfy one of the following two conditions:

   1. {ξₖ} is uniformly bounded and ηₖ=o((1)/(log k)).
   2. {ξₖ} has bounded variance, i.e. 𝔼[‖ ξₖ₊₁ ‖ ²|Fₖ]≤σ²<∞, and Σₖ₌₀^∞ ηₖ²<∞.

4. limₖ→∞ δₖ=0 and limₖ→∞ νₖ=0.

5. There exists a positive τ such that limₖ→∞(θₖ)/(ηₖ)=τ.

We make some remarks on Assumption 4.3. Assumption 4.3(1) and (3) are directly adapted from Assumption 3.1(1) and (4), respectively. Assumption 4.3(4) characterizes how dₖ approximates D_f(xₖ) and quantifies the inexactness in solving the subproblem. We note that Assumption 4.3(2) can be ensured by the global stability when the noise sequence {ξₖ} corresponds to randomly reshuffling noises, as demonstrated by Theorem 3.3. Assumption 4.3(5) assumes that the stepsizes {ηₖ} and momentum parameters {θₖ} in the framework (MSBG) are single-timescale, i.e. they decrease at the same rate. In practice, τ can be preset.

Consider the following differential inclusion,

$$\frac{d}{dt}\left[\begin{array}{c}\nabla\varphi(x(t))\\ m(t)\end{array}\right]\in-\left[\begin{array}{c}m(t)\\ \tau(\nabla^{2}\varphi(x(t)))^{-1}(m(t)-D_{f}(x(t)))\end{array}\right],\text{for almost all}t\geq 0\text{.}$$

Define

$$\phi(x,m):=\varphi(x)+\frac{1}{2}\Vert m\Vert ^{2},\;\quad\mathcal{H}(x,m):=\left[\begin{array}{c}m\\ \tau(\nabla^{2}\varphi(x))^{-1}(m-D_{f}(x))\end{array}\right].$$

Then, (17) can be reformulated in the form of (6) as:

$$\frac{d}{dt}\nabla\phi(x(t),m(t))\in-\mathcal{H}(x(t),m(t)),\text{for almost all}t\geq 0\text{.}$$

(18)

The stable set of (18) is given by H⁻¹(0)={(x,m): m=0,0∈ D_f(x)}. Based on the differential inclusion (18), (MSBG) can be reformulated in the form of (8) as:

$$\left\lbrace \nabla\varphi(x_{k+1})=\nabla\varphi(x_{k})-\eta_{k}(m_{k}+d^{e}_{k})\quad m_{k+1}=m_{k}-\eta_{k}\cdot\frac{\theta_{k}}{\eta_{k}}(\nabla^{2}\varphi(x_{k}))^{-1}(m_{k}-d_{k}), \right.$$

(19)

where ‖ dᵉₖ ‖≤νₖ.

Lemma 6

Suppose Assumptions 4.3 and 4.2 hold. Let {(xₖ,mₖ)} be the sequence generated by (MSBG), dₓ,ₖ:=mₖ+dᵉₖ, and dₘ,ₖ:=(θₖ)/(ηₖ)(∇² φ(xₖ))⁻¹(mₖ-dₖ), where ‖ dᵉₖ ‖≤νₖ. For any increasing sequence {kⱼ} such that (x_kⱼ,m_kⱼ) converges to (x̅,m̅), it holds that

$$\lim_{N\rightarrow\infty}\mathrm{dist}\left(\frac{1}{N}\sum_{j=1}^{N}(d_{x,k_{j}},d_{m,k_{j}}),\mathcal{H}(\bar{x},\bar{m})\right)=0.$$

Proof

By Assumption 4.3, it holds that

$$\lim_{j\rightarrow\infty}\mathrm{dist}\left(m_{k_{j}}+d^{e}_{k_{j}},\bar{m}\right)\leq\lim_{j\rightarrow\infty}\mathrm{dist}(m_{k_{j}},\bar{m})+\nu_{k_{j}}=0.$$

By Assumption 4.3.1, we have (∇² φ(x))⁻¹ is continuous, and hence

$$\lim_{j\rightarrow\infty}\mathrm{dist}\left((\nabla^{2}\varphi(x_{k_{j}}))^{-1}(m_{k_{j}}-d_{k_{j}}),(\nabla^{2}\varphi(\bar{x}))^{-1}(\bar{m}-D_{f}(\bar{x}))\right)=0.$$

Since D_f(x̅) is a compact set, and limₖ→∞(θₖ)/(ηₖ)=τ, it holds that

$$\lim_{j\rightarrow\infty}\mathrm{dist}\left(\frac{\theta_{k_{j}}}{\eta_{k_{j}}}(\nabla^{2}\varphi(x_{k_{j}}))^{-1}(m_{k_{j}}-d_{k_{j}}),\tau(\nabla^{2}\varphi(\bar{x}))^{-1}(\bar{m}-D_{f}(\bar{x}))\right)=0.$$

By the fact that D_f(x̅) is a convex set and Jensen’s inequality, we have that

$$\lim_{N\rightarrow\infty}\mathrm{dist}\left(\frac{1}{N}\sum_{j=1}^{N}(d_{x,k_{j}},d_{m,k_{j}}),\mathcal{H}(\bar{m},\bar{x})\right)\leq\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{j=1}^{N}\mathrm{dist}\left((d_{x,k_{j}},d_{m,k_{j}}),\mathcal{H}(\bar{m},\bar{x})\right)=0.$$

This completes the proof.

Proposition 2

Suppose Assumptions 4.3 and 4.2 hold. Then h(x,m)=f(x)+(1)/(2τ)‖ m‖ ² is a Lyapunov function for (17) with the stable set B:={(x,m)∈ℝⁿ ×ℝⁿ: m=0,0∈ D_f(x)}.

Proof

Consider any trajectory (x(t),m(t)) for the differential inclusion (17) with (x(0),m(0))∉B. There exists measurable d_f(s)∈ D_f(x(s)), such that for almost all s≥0,

$$\frac{d}{ds}h(x(s),m(s))\qquad {} = \left\langle D_{f}(x(s)),\dot{x}(s)\right\rangle+\left\langle\frac{m(s)}{\tau},\dot{m}(s)\right\rangle\qquad \ni -\left\langle d_{f}(s),(\nabla^{2}\varphi(x(s)))^{-1}m(s)\right\rangle-\left\langle m(s),(\nabla^{2}\varphi(x))^{-1}(m(s)-d_{f}(s))\right\rangle\qquad {} = -\left\langle m(s),(\nabla^{2}\varphi(x(s)))^{-1}m(s)\right\rangle.$$

Thus, for any t≥0, h(x(t),m(t))≤h(x(0),m(0)). For any (x(0),m(0))∉B, either m(0)≠0 or m(0)=0 and 0∉D_f(x(0)). If m(0)≠0, then the continuity of m(·) ensures the existence of T>0 and α>0 where ‖ m(s)‖≥α for s∈[0,T]. Thus we have

$$h(x(T),m(T))-h(x(0),m(0))\leq-\int_{0}^{T}\left\langle m(s),(\nabla^{2}\varphi(x(s)))^{-1}m(s)\right\rangle ds\lt 0.$$

Now consider the case m(0)=0 and 0∉D_f(x(0)). By the outer semicontinuity of D_f and Assumption 4.3.1, there exists T̃>0, such that for almost all t∈[0,T̃], it holds that 0∉(∇² φ(x(t)))⁻¹ D_f(x(t)). Now suppose for all t≥0, h(x(t),m(t))=h(x(0),m(0)), then we have m(s)=0 for almost all s≥0. Since m is continuous, it holds that m≡0. Note that for almost any t≥0, ṁ(t)∈-τ(∇² φ(x(t)))⁻¹(m(t)-D_f(x(t))), thus 0∈(∇² φ(x(t)))⁻¹ D_f(x(t)) holds for almost all t≥0, which leads to a contradiction. Therefore, for both cases, there exists T>0 such that h(x(T),m(T))<h(x(0),m(0)). This completes the proof.

By Lemma 6, Proposition 2, and Theorem 3.2, we have the following convergence results for (MSBG).

Theorem 4.2

Suppose Assumptions 4.3 and 4.2 hold. Then almost surely, any cluster point of {xₖ} generated by (MSBG) is a D_f-stationary point of f, limₖ→∞ mₖ=0, and the function values {f(xₖ)} converge.

Proof

Theorem 3.2 implies that any cluster point of {(xₖ,mₖ)} lies in {(x,m):0∈ D_f(x),m=0}, and {f(xₖ)+(1)/(2τ)‖ mₖ ‖ ²} converges. For any convergent subsequence x_k_j→x, since {mₖ} is bounded, then there exist subsequence {m_k_jᵢ} such that m_k_j_i→m. Therefore, (x_k_jᵢ,m_k_jᵢ)→(x̅,m̅). Then, it holds that 0∈ D_f(x̅). Similarly, we can prove that for any convergent subsequence {m_kⱼ} such that m_k_j→m, we have that m̅=0. Therefore, limₖ→∞ mₖ=0. Then, limₖ→∞ f(xₖ)=limₖ→∞ f(xₖ)+(1)/(2τ)‖ mₖ ‖ ². This completes the proof.

4.3 Stochastic Bregman proximal subgradient method

In this section, we consider solving the following constrained composite optimization problem:

$$\min_{x\in\mathcal{X}}\;h(x):=f(x)+R(x),$$

(20)

where X is a closed subset of ℝⁿ, and R is a locally Lipschitz function with an efficiently computable conservative field. In many applications, R serves as the regularization function, which is usually Clarke regular, and ∂ R is efficient to compute. We consider applying the follow Bregman proximal subgradient method to solve (20),

x_k+1 argmin _x∈ X ⟨ g_k,x-x_k⟩+1_k D _ϕ(x,x_k)+R(x) , s.t. ⟨ g_k,x_k+1-x_k⟩+1_k D ϕ(x_k+1,x_k)+R(x_k+1)≤ R(x_k), and dist (0,g_k+1_k(ϕ(x_k+1)-ϕ(x_k))+D_R(x_k+1)+N^L X (x_k+1))≤ν_k. .

(SBPG)

where gₖ=d_[f,k]+ξₖ, d_[f,k] ∈ Dᵟₖ_f(xₖ). We can reformulate (SBPG) in the form of (8) as follows:

$$\nabla\phi(x_{k+1})=\nabla\phi(x_{k})-\eta_{k}(d_{f,k}+d_{R,k}+d_{\mathcal{X},k}+d_{e,k}+\xi_{k}),$$

where d_[f,k] ∈ Dᵟₖ_f(xₖ), d_[R,k] ∈ D_R(xₖ₊₁), d_X,k∈ N^L_X(xₖ₊₁), and ‖ dₑ,ₖ ‖≤νₖ. When X=ℝⁿ, with δₖ=0 and νₖ=0, it follows that xₖ₊₁=T^ϕ_ηₖ,R(xₖ,gₖ). Let dₖ:=d_[f,k]+d_[R,k]+d_X,k+dₑ,ₖ. We can easily verify that there exists {δ̃ₖ} such that limₖ→∞δ̃ₖ=0, and dₖ ∈Hᵟ̃ₖ(xₖ), where H:=D_f+D_R+N_X. This leads to a differential inclusion for the proximal updates given by

$$\frac{d}{dt}\nabla\phi(x(t))\in-\mathcal{H}(x(t)),\;\text{where}\mathcal{H}=D_{f}+D_{R}+N_{\mathcal{X}}.$$

(21)

The momentum technique can also be integrated into (SBPG), as illustrated in Section 4.2. For the sake of readability, we omit this extension. Note that neither N_X nor N^L_X is locally bounded, thus the results that rely on local boundedness assumption such as those presented in xiao2023convergence ; ding2023adam cannot be directly applied. We make the following assumptions on (SBPG).

Assumption 4.4

1. ϕ is a supercoecive Legendre kernel function over ℝⁿ. Moreover, ϕ is locally strongly convex and ∇ϕ is locally Lipschitz continuous. Additionally, for any absolutely continuous mapping z(·)∈C(ℝ₊,ℝⁿ), ∇² ϕ(z(s)) is positive definite for almost all s≥0.
2. The sequences {xₖ}, {∇ϕ(xₖ)}, {d_[f,k]} and {d_[R,k]} are uniformly bounded almost surely.
3. {ξₖ} is a uniformly bounded martingale difference noise and supₖ≥0 ‖ ξₖ ‖<∞. The stepsize sequence {ηₖ} satisfies Σₖ₌₀^∞ ηₖ=∞ and ηₖ=o((1)/(log k)).
4. limₖ→∞ δₖ=0, and limₖ→∞ νₖ=0.
5. For any η>0, η R+ϕ is supercoecive.

Now, we make some remarks on Assumption 4.4. Assumption 4.4(1) and (3) are directly adapted from Assumption 3.1(1) and (4), respectively. Assumption 4.1(4) characterizes how dₖ approximates D_f(xₖ) and quantifies the inexactness in solving the subproblem. Assumption 4.1(2) can be ensured by the global stability when the noise sequence {ξₖ} corresponds to randomly reshuffling noises, as demonstrated by Theorem 3.3. Assumption 4.4(5) ensures the well-posedness of Bregman proximal mapping as demonstrated in Lemma 3.

To ensure the convergence of (SBPG), we make the following assumptions on f and kernel ϕ, which is essentially a specific case of Assumption 3.2.

Assumption 4.5

1. h is lower bounded. Moreover, the locally Lipschitz continuous functions f and R are potential functions that admit convex-valued D_f and D_R as their conservative fields, respectively. X admits the chain rule as described in Definition 6.
2. The critical value set {h(x):0∈ D_f(x)+D_R(x)+N_X(x)} has empty interior in ℝ.

Lemma 7

Suppose Assumptions 4.4 and 4.5 hold. Let dₖ:=d_[f,k]+d_[R,k]+d_X,k+dₑ,ₖ. For any increasing sequence {kⱼ} such that {x_kⱼ} converges to x̅, it holds that

$$\lim_{N\rightarrow\infty}\mathrm{dist}\left(\frac{1}{N}\sum_{j=1}^{N}d_{k_{j}},\mathcal{H}(\bar{x})\right)=0,$$

where H:=D_f+D_R+N_X.

Proof

Given the condition

$$R(x_{k+1})+\left\langle d_{f,k}+\xi_{k},x_{k+1}-x_{k}\right\rangle+\frac{1}{\eta_{k}}\mathcal{D}_{\phi}(x_{k+1},x_{k})\leq R(x_{k}),$$

we derive that

$$\frac{\mathcal{D}_{\phi}(x_{k+1},x_{k})}{\left\Vert x_{k+1}-x_{k}\right\Vert }\leq\eta_{k}\frac{|R(x_{k+1})-R(x_{k})|}{\left\Vert x_{k+1}-x_{k}\right\Vert }+\eta_{k}(\left\Vert d_{f,k}+\xi_{k}\right\Vert ).$$

Assumption 4.4 ensures that supₖ≥0(1)/(ηₖ)‖ ∇ϕ(xₖ₊₁)-∇ϕ(xₖ)‖<∞. Moreover, we have that for some ‖ dₑ,ₖ ‖≤νₖ,

$$\frac{\nabla\phi(x_{k+1})-\nabla\phi(x_{k})}{\eta_{k}}=-\left(d_{f,k}+d_{R,k}+d_{\mathcal{X},k}+d_{e,k}+\xi_{k}\right).$$

Note that the left hand side is uniformly bounded, {d_[f,k]}, {d_[R,k]}, {dₑ,ₖ} and {ξₖ} are all uniformly bounded, therefore, it holds that {d_X,k} is also uniformly bounded. For any {x_kⱼ} converging to x̅, by the outer semicontinuity of D_f, D_R and N_X^L, along with (RockWets98, , Theorem 5.7, Proposition 5.12), it holds that

$$\lim_{j\rightarrow\infty}\mathrm{dist}\left(d_{x,k_{j}},D_{f}(\bar{x})\right)=0,\;\lim_{j\rightarrow\infty}\mathrm{dist}\left(d_{R,k_{j}},D_{R}(\bar{x})\right)=0,\;\text{and}\lim_{j\rightarrow\infty}\mathrm{dist}\left(d_{\mathcal{X},k_{j}},N^{L}_{\mathcal{X}}(\bar{x})\right)=0.$$

Note that N^L_X⊂ N_X, therefore, limⱼ→∞ dist(d_kⱼ,H(x̅))=0. By Jensen’s inequality, we prove this lemma.

Proposition 3

Suppose Assumptions 4.4 and 4.5 hold. Then the function h in (20) is a Lyapunov function for the differential inclusion (21) with stable set {x∈ℝⁿ:0∈ D_f(x)+D_R(x)+N_X(x)}.

Proof

Consider any trajectory z(t) of (21) with 0∉H(z(0)). By the chain rule, it holds that for almost all t≥0,

$$f(z(t))^{\prime}=\left\langle D_{f}(z(t)),\dot{z}(t)\right\rangle,\;R(z(t))^{\prime}=\left\langle D_{R}(z(t)),\dot{z}(t)\right\rangle,\;0=\left\langle N_{\mathcal{X}}(z(t)),\dot{z}(t)\right\rangle.$$

Then, we have that h(z(t))′=⟨H(z(t)),ż(t)⟩ for almost all t≥0. Note that for almost all t≥0, it holds that ∇² ϕ(z(t))ż(t)∈-H(z(t)). Then, we have

$$h(z(t))-h(z(0))=\int_{0}^{t}\left\langle\mathcal{H}(z(s)),\dot{z}(s)\right\rangle ds=-\int_{0}^{t}\left\langle\nabla^{2}\phi(z(s))\dot{z}(s),\dot{z}(s)\right\rangle ds\qquad {} \leq -\int_{0}^{t}\lambda_{\min}(\nabla^{2}\phi(z(s)))\left\Vert \dot{z}(s)\right\Vert ^{2}ds.$$

If there exists no t>0, such that h(z(t))<h(z(0)). Then, we have that ż(t)=0 for almost all t≥0. Thus, z(t)≡z(0). Therefore, we have that 0=(d)/(dt)∇ϕ(z(s))∈-H(z(0)), which is contradictory to the fact that 0∉H(z(0)). This completes the proof.

By Lemma 7, Proposition 3 and Theorem 3.2, we can directly derive the following convergence results.

Theorem 4.3

Suppose Assumptions 4.4 and 4.5 hold. Then almost surely, any cluster point of {xₖ} generated by (SBPG) is a H-stationary point and the function values {f(xₖ)+R(xₖ)} converge.

5 Numerical experiments

In this section, we conduct preliminary numerical experiments to illustrate the performance of our proposed methods, focusing on training nonsmooth neural networks for image classification and language modeling tasks. These experiments are performed using an NVIDIA RTX 3090 GPU, and implemented in Python 3.9 with PyTorch version 1.12.0.

Our experiments employ a polynomial kernel-based stochastic Bregman subgradient method to train nonsmooth neural networks. Specifically, we use a blockwise polynomial kernel function ϕ(x)=Σᵢ₌₁^L pᵢ(‖ xᵢ ‖), where x=(x₁,...,x_L) represents the concatenation of all layers’ parameters in a neural network with L layers, and each pᵢ is a univariate polynomial of degree at least 2. When pᵢ(λ)=1/2λ², this approach becomes equivalent to SGD. The polynomial pᵢ(λ)=1/2λ²+(σ)/(r)λʳ with r≥4, as discussed in the prior work ding2023nonconvex , is applied in our numerical experiments. In this case, the update scheme is defined as follows

$$x_{k+1} =\nabla\phi^{\ast }(\nabla\phi(x_{k})-\eta_{k}m_{k})\qquad m_{k+1} =m_{k}-\theta_{k}(\nabla^{2}\phi(x_{k}))^{-1}(m_{k}-g_{k}),$$

where the calculation of xₖ₊₁ involves solving a nonlinear equation. Given the allowance of inexact solutions for the subproblems of MSBG, as mentioned earlier, we consider an alternative preconditioned update scheme,

$$x_{k+1} =x_{k}-\eta_{k}(\nabla^{2}\phi(x_{k}))^{-1}m_{k}\qquad m_{k+1} =m_{k}-\theta_{k}(\nabla^{2}\phi(x_{k}))^{-1}(m_{k}-g_{k}),$$

which avoids solving a nonlinear equation in the first step of (MSBG). By the Sherman-Morrison formula, we have

$$(\nabla|_{x_{i}}^{2}p_{i}(\Vert x_{i}\Vert ))^{-1}=\frac{1}{\sigma\left\Vert x_{i}\right\Vert ^{r-2}}I-\frac{\sigma(r-2)\left\Vert x_{i}\right\Vert ^{r-4}}{(1+\sigma\left\Vert x_{i}\right\Vert ^{r-2})^{2}+\sigma(r-2)(1+\sigma\left\Vert x_{i}\right\Vert ^{r-2})\left\Vert x_{i}\right\Vert ^{2}}x_{i}x_{i}^{T},$$

and ∇² ϕ(x)=diag((∇|_x₁² p₁(‖ x₁ ‖))⁻¹,...,(∇|_x_L² p_L(‖ x_L ‖))⁻¹) is block diagonal. We employ this kernel function and use the notation MSBGK/iMSBGK to denote MSBG/iMSBG with the polynomial degree parameter r set to K. Our experiments focus on two main applications: training Convolutional Neural Networks (CNNs) for image classification and Long Short-Term Memory (LSTM) hochreiter1997long networks for language modeling. Specifically, our image classification experiments include training Resnet14 and ResNet34 he2016deep on CIFAR-10 and CIFAR-100 datasets krizhevsky2009learning . Our language modeling experiments focus on 1-layer, 2-layer, and 3-layer LSTM networks applied to the Penn Treebank dataset marcus1993building .

CNNs on image classification

For the CNN experiments, we set the stepsize ηₛ for each epoch s as ηₛ=(η₀)/(1+(log(s+1))¹.¹), where η₀ is the initial stepsize. The momentum parameters are all set to θₛ=(0.1)/(1+log(s+1)¹.¹). For MSBG4 and iMSBG4, we choose σ=0.01, and for iMSBG6, σ=0.0001. We search the initial stepsize η₀ among the grid {0.001,0.01,0.1,1.0} and select the value that achieves the highest test accuracy. The results are shown in Figure 1 and Figure 2. We can observe that by selecting a proper kernel function, our Bregman subgradient methods can outperform SGD in terms of test accuracy. Moreover, we can see that MSBG4 and iMSBG4 have similar performance, although iMSBG4 solves the subproblem inexactly.

Additionally, we evaluate the robustness of the selection of initial stepsize η₀, as demonstrated in Figure 3. We can see that the peak test accuracies of all methods are similar, yet our Bregman subgradient methods demonstrate a wider effective initial step size range, indicating a reduced sensitivity to the choice of initial step size – a benefit attributable to the kernel function.

Figure 1: Resnet 14 on CIFAR10.

Figure 2: Resnet 34 on CIFAR100.

Figure 3: Robustness for initial stepsize. Figure (a) reports the test accuracy for 40 epochs. Figure (b) reports the test accuracy for 30 epochs.

LSTMs on language modeling

For the LSTM experiments, we initially set the stepsize as a constant. The stepsize is then decreased to 0.1 times its previous value at both the 150th and 300th epochs. After 300 epochs, we set ηₛ=(0.01η₀)/(1+log(s-300)¹.¹), with s representing the epoch number. Here η₀ is the initial stepsize. Within the s-th epoch, ηₖ takes the constant value ηₛ. Similarly, the momentum parameters are all set to θₛ=(0.1)/(1+log(s+1)¹.¹). For MSBG4 and iMSBG4, we set σ=10⁻⁶. We search η₀ among the grid {1,10,20,40,80,100} and report the results based on achieving the highest test accuracy. The results are shown in Figures 4, 5, and 6. We can observe that selecting an appropriate kernel function enables our Bregman subgradient methods to achieve superior test accuracy compared to SGD.

We also compare the one-epoch runtime for all considered methods over all tasks. We can observe in Table 1 that the proposed inexact Bregman subgradient methods are nearly as efficient as SGD, largely because iMSBG circumvents the need to solve nonlinear equations in computing the Bregman proximal mapping.

Figure 4: 1-layer LSTM on Penn Treebank dataset.

Figure 5: 2-layer LSTM on Penn Treebank dataset.

Figure 6: 3-layer LSTM on Penn Treebank dataset.

Table 1: Computation Time of Each Epoch (in seconds). Task 1 is training ResNet14 on CIFAR10. Task 2 is training ResNet18 on CIFAR100. Task 3,4,5 are training one-layer, two-layer, three-layer LSTM on n Penn Treebank dataset, respectively.

6 Conclusions

This paper explores Bregman subgradient methods for solving nonsmooth nonconvex optimization problems, particularly focusing on path-differentiable functions. We introduce a comprehensive stochastic Bregman framework that accommodates inexact evaluations of the abstract set-valued mapping. Employing a differential inclusion strategy and linear interpolation of dual sequences, we establish convergence results for our stochastic Bregman-type methods. This ensures that the discrete sequence subsequentially converges to the stable set of the differential inclusion, as well as the convergence of the Lyapunov function values. For applications, we demonstrate that stochastic Bregman subgradient methods, even with subproblems being solved inexactly, fit within our general framework, and we establish their convergence properties. Moreover, we integrate a momentum technique into the stochastic Bregman subgradient methods. Additionally, we extend our methodology to a proximal variant of the stochastic Bregman subgradient methods for solving constrained composite optimization problems and establish its convergence results. Finally, we conduct numerical experiments on training nonsmooth neural networks to evaluate the performance of our proposed stochastic Bregman subgradient methods. Our experimental results validate the practical benefits and effectiveness of our approaches in deep learning.

Funding The research of Kim-Chuan Toh is supported by the Ministry of Education, Singapore, under its Academic Research Fund Tier 2 grant (MOE-T2EP20224-0029).

Data Availibility The datasets analyzed during the current study are publicly available in the following resources:

• The CIFAR-10 and CIFAR-100 datasets, introduced by krizhevsky2009learning , can be accessed online at https://www.cs.toronto.edu/~kriz/cifar.html.
• The Penn Treebank dataset, described by marcus1993building , is available at https://catalog.ldc.upenn.edu/LDC99T42.

Appendix Appendix A: Proof of Theorem 3.1

To establish under Assumption 3.1 that, the discrete sequence generated by the general Bregman-type method (8) tracks a trajectory of the differential inclusion (6), we introduce a piecewise constant mapping d(·) defined by d(s)=dₖ for any s∈[λ_η(k),λ_η(k+1)). We also define a time-shifted solution xₜ(·) to the following ordinary differential equation:

$$\frac{d}{ds}\nabla\phi({x_{t}(s)})=-d(s)\;for\;all\;s\geq t,\;with initial condition\;\nabla\phi(x_{t}(t))=\nabla\phi(x(t)).$$

The following lemma suggests that the interpolated process x(t) defined by (9) asymptotically approximates this time-shifted solution. We remind the reader to note the difference between xₜ(·) and xᵗ(·).

Lemma A.8

Suppose Assumption 3.1 holds, then for any T>0, it holds that

$$\lim_{t\rightarrow\infty}\sup_{s\in[t,t+T]}\Vert \nabla\phi(x(s))-\nabla\phi(x_{t}(s))\Vert =0.$$

(22)

Proof

Fix an arbitrary s∈[t,t+T], let τₜ=Λ_η(t), τₛ=Λ_η(s). By the definition of xₜ(·), we have

∇ϕ(xₜ(s))=∇ϕ(x(t))-∫ₜ^λ_η(τₜ)d(u)du-∫_λ_η(τₜ)^λ_η(τₛ)d(u)du-∫_λ_η(τₛ)ˢ d(u)du=∇ϕ(x(λ_η(τₜ)))-Σ_i=τₜ^τₛ-1ηᵢ dᵢ+(∇ϕ(x(t))-∫ₜ^λ_η(τₜ)d(u)du-∇ϕ(x(λ_η(τₜ))))-∫_λ_η(τₛ)ˢ d(u)du=∇ϕ(x(λ_η(τₛ)))+Σ_i=τₜ^τₛ-1ηᵢ ξᵢ+(∇ϕ(x(t))-∫ₜ^λ_η(τₜ)d(u)du-∇ϕ(x(λ_η(τₜ))))-∫_λ_η(τₛ)ˢ d(u)du=∇ϕ(x(s))+Σ_i=τₜ^τₛ-1ηᵢ ξᵢ+(∇ϕ(x(t))-∫ₜ^λ_η(τₜ)d(u)du-∇ϕ(x(λ_η(τₜ))))

$${}+\left(\nabla\phi(x(\lambda_{\eta}(\tau_{s})))-\int_{\lambda_{\eta}(\tau_{s})}^{s}d(u)du-\nabla\phi(x(s))\right).$$

Note that

$${} \left\Vert \nabla\phi(x(t))-\int_{t}^{\lambda_{\eta}(\tau_{t})}d(u)du-\nabla\phi(x(\lambda_{\eta}(\tau_{t})))\right\Vert\qquad {} \leq {\Vert \nabla\phi(x(t))-\nabla\phi(x(\lambda_{\eta}(\tau_{t})))\Vert }+\int_{t}^{\lambda_{\eta}(\tau_{t})}\Vert d(u)\Vert du\qquad {} \leq {\Vert \nabla\phi(x(\lambda_{\eta}(\tau_{t}+1)))-\nabla\phi(x(\lambda_{\eta}(\tau_{t})))\Vert }+\int_{t}^{\lambda_{\eta}(\tau_{t})}\Vert d(u)\Vert du\qquad {} \leq \eta_{\tau_{t}}(\Vert \xi_{\tau_{t}}\Vert +2\Vert d_{\tau_{t}}\Vert ),$$

and similarly ‖ ∇ϕ(x(λ_η(τₛ)))-∫_λ_η(τₛ)ˢ d(u)du-∇ϕ(x(s))‖≤η_τₛ(‖ ξ_τₛ‖ +2‖ d_τₛ‖). By Assumption 3.1, we have that lim_t→∞_τ_t(‖ _τ_t‖ +2‖ d_τ_t‖ )=0, lim_s→∞_τ_s(‖ _τ_s‖ +2‖ d_τ_s‖ )=0, and lim_t→∞_s∈[t,t+T]Σ_i=τ_t^τ_s-1_i_i=0. Therefore, it holds that lim_t→∞_s∈[t,t+T]‖ ϕ(x(s))-ϕ(x_t(s))‖ =0, which completes the proof.

Proof of Theorem 3.1 . By the definition of xₜ(·), it follows that ∇ϕ(xₜ(s))=∇ϕ(x(t))-∫ₜˢ y(u)du, for all s≥t. By the boundedness of y(s), Arzela-Ascoli’s theorem rudin1953principles ensures that {∇ϕ(xₜ(t+·))} _t∈ℝ₊ is relatively compact in C(ℝ₊,ℝⁿ). For any subset {τₖ} ⊂ℝ₊, we consider the sequence {∇ϕ(x^τₖ(·))} _k∈ℕ₊. There are two cases to consider. Case (i): the sequence {τₖ} has a cluster point t. Without loss of generality, assume that limₖ→∞ τₖ=t. By the definition of x(·) in (9) and uniform boundedness assumption, it holds that ∇ϕ(x(·)) is Lipschitz continuous. Thus, ∇ϕ(x^τₖ(·)) converges to ∇ϕ(xᵗ(·)) in C(ℝ₊,ℝⁿ). Case (ii): limₖ→∞ τₖ=∞. Suppose that ∇ϕ(x^τₖ(·)) does not has any cluster point in C(ℝ₊,ℝⁿ). Since {∇ϕ(x_τₖ(τₖ+·))} is relatively compact in C(ℝ₊,ℝⁿ), without loss of generality, we assume that limₖ→∞ ∇ϕ(x_τₖ(τₖ+·))=y̅(·). Then, for any compact set C⊂ℝ₊, it follows from Lemma A.8 that

$$\lim_{k\rightarrow\infty,s\in C}\left\Vert \nabla\phi(x^{\tau_{k}}(s))-\bar{y}(s)\right\Vert\qquad {} \leq \lim_{k\rightarrow\infty,s\in C}\left\Vert \nabla\phi(x^{\tau_{k}}(s))-\nabla\phi(x_{\tau_{k}}(\tau_{k}+s))\right\Vert +\lim_{k\rightarrow\infty,s\in C}\left\Vert \nabla\phi(x_{\tau_{k}}(\tau_{k}+s))-\bar{y}(s)\right\Vert =0,$$

which contradicts that ∇ϕ(x^τₖ(·)) does not has any cluster point in C(ℝ₊,ℝⁿ). Thus, in both cases, {∇ϕ(x^τₖ(·))} is relatively compact in C(ℝ₊,ℝⁿ). Because {τₖ} is an arbitrary subset in ℝ₊, we have that {∇ϕ(xᵗ(·))} is relatively compact in C(ℝ₊,ℝⁿ).

Next, we aim to construct a trajectory of the differential inclusion. Define the shifts dᵗ(·)=d(t+ ·). Consider {τₖ} satisfying τₖ →∞, and fix T>0. Without loss of generality, we assume that ∇ϕ(x^τₖ(·)) converges to y̅(·) in C(ℝ₊,ℝⁿ), otherwise, we choose its convergent subsequence. The set Y_T:={d^τₖ(s),s∈[0,T]}ₖ∈ℕ ⊂ L²([0,T]) is bounded. Therefore, it follows from the Banach-Alaoglu theorem rudin1991functional that Y_T is weakly sequentially compact, i.e. there exists a subsequence {τ_kⱼ} and y̅(·)∈ L²([0,T]) such that d^τ_k_j(·)→d(·) weakly in L²([0,T]). On the other hand, by Lemma A.8, we have ∇ϕ(x_τ_kⱼ(τ_kⱼ+·)) converges to y̅(·) in C(ℝ₊,ℝⁿ). For any τ∈[0,T], by the definition of xₜ(·), we have

$$\nabla\phi(x_{t}(t+\tau))=\nabla\phi(x_{t}(t))-\int_{0}^{\tau}d^{t}(s)ds.$$

Setting t=τ_kⱼ and taking the limit as kⱼ →∞, we deduce that

$$\bar{y}(\tau)=\bar{y}(0)-\int_{0}^{\tau}\bar{d}(s)ds.$$

Let x̅(·)=∇ϕ∗(y̅(·)), since T>0 is arbitrary, we get (12).

The remaining step is to verify that d̅(s)∈H(x̅(s)) for almost all s≥0. We again fix an arbitrary T>0. Given that Y_T ⊂ L²([0,T]) is bounded, the Banach-Saks theorem rudin1991functional implies that for τₖ (choosing a subsequence if necessary), (1)/(N)Σₖ₌₁^N d^τₖ(s) strongly converges to d̅(s) in L²([0,T]) . By the definition of d(·), we have d^τₖ(s)=d_Λ_η(τₖ+s). Now for any s∈[0,T], we have

$$\Vert \nabla\phi(x(\lambda_{\eta}(\Lambda_{\eta}(\tau_{k}+s))))-\nabla\phi(\bar{x}(s))\Vert\qquad {} \leq \Vert \nabla\phi(x(\lambda_{\eta}(\Lambda_{\eta}(\tau_{k}+s))))-\nabla\phi(x(\tau_{k}+s))\Vert +\Vert \nabla\phi(x^{\tau_{k}}(s))-\nabla\phi(\bar{x}(s))\Vert\qquad {} \leq \Vert \nabla\phi(x(\lambda_{\eta}(\Lambda_{\eta}(\tau_{k}+s)+1)))-\nabla\phi(x(\lambda_{\eta}(\Lambda_{\eta}(\tau_{k}+s))))\Vert +\Vert \nabla\phi(x^{\tau_{k}}(s))-\nabla\phi(\bar{x}(s))\Vert\qquad {} \leq \eta_{\Lambda_{\eta}(\tau_{k}+s)}(\Vert \xi_{\Lambda_{\eta}(\tau_{k}+s)}\Vert +\Vert d_{\Lambda_{\eta}(\tau_{k}+s)}\Vert )+\Vert \nabla\phi(x^{\tau_{k}}(s))-\nabla\phi(\bar{x}(s))\Vert ,$$

which converges to zero as k→∞. By the continuity of ∇ϕ∗, we have that x(λ_η(Λ_η(τₖ+s))) converges to x̅(s) in C(ℝ₊,ℝⁿ). By Assumption 3.1, for almost any s∈[0,T], we have

dist (d(s), H (x(s))) ≤‖ (1)/(N)Σ_k=1^Nd^τ_k(s)-d(s)‖ + dist ((1)/(N)Σ_k=1^Nd^τ_k(s), H (x(s))) = ‖ (1)/(N)Σ_k=1^Nd^τ_k(s)-d(s)‖ + dist ((1)/(N)Σ_k=1^Nd__(λ_(τ_k+s)), H (x(s)))→ 0.

Since T is arbitrary and H(x̅(s)) is a closed set, we conclude that d̅(s)∈H(x̅(s)) for almost all s≥0. This completes the proof.∎

Appendix Appendix B: Proof of Theorem 3.2

Lemma B.9

Suppose Assumption 3.1 is satisfied. Then, it holds that

$$\liminf_{t\rightarrow\infty}\Psi(x(t))=\liminf_{k\rightarrow\infty}\Psi(x_{k}),\quad\limsup_{t\rightarrow\infty}\Psi(x(t))=\limsup_{k\rightarrow\infty}\Psi(x_{k}).$$

(23)

Proof

For simplicity, we only prove the case for $\liminf$, the argument for $\limsup$ follows similarly. By Assumption 3.1, we have that

$$\lim_{k\rightarrow\infty}\Vert \nabla\phi(x_{k+1})-\nabla\phi(x_{k})\Vert =0.$$

(24)

Let τᵢ →∞ be an arbitrary sequence with x(τᵢ)→ x∗. By the definition of λ_η and Λ_η, it follows that

$$\Vert \nabla\phi(x_{\Lambda_{\eta}(\tau_{i})})-\nabla\phi(x^{\ast })\Vert \leq \Vert \nabla\phi(x_{\Lambda_{\eta}(\tau_{i})})-\nabla\phi(x(\tau_{i}))\Vert +\Vert \nabla\phi(x(\tau_{i}))-\nabla\phi(x^{\ast })\Vert\qquad {} \leq \Vert \nabla\phi(x_{\Lambda_{\eta}(\tau_{i})})-\nabla\phi(x_{\Lambda_{\eta}(\tau_{i})+1})\Vert +\Vert \nabla\phi(x(\tau_{i}))-\nabla\phi(x^{\ast })\Vert .$$

The right-hand side converges to zero. By Remark 1.1, we have that ϕ^ is continuous, so lim_i→∞(x__(τ_i))=lim_i→∞ϕ^ ϕ(x__(τ_i))=ϕ^ ϕ(x^ )=(x^ ). By choosing τᵢ →∞ as the sequence realizing _t→∞(x(t)), and assuming without loss of generality that x(τᵢ)→ x∗, we get

$$\liminf_{k\rightarrow\infty}\Psi(x_{k})\leq\lim_{{i\rightarrow\infty}}\Psi(x_{\Lambda_{\eta}(\tau_{i})})=\Psi(x^{\ast })=\liminf_{t\rightarrow\infty}\Psi(x(t)).$$

This completes the proof.

The following lemma demonstrates that the function value converges along the interpolated process defined in (9). The non-escape argument in the proof is adapted from those of (davis2020stochastic, , Proposition 3.5) and (duchi2018stochastic, , Theorem 3.20), with particular attention paid to the dual map ∇ϕ and its inverse ∇ϕ∗.

Lemma B.10

Suppose Assumption 3.1 and 3.2 hold, then function value Ψ(x(t)) converges as t→∞.

Proof

Assuming $\liminf_{t\rightarrow\infty}\Psi(x(t))=0$, we define the level set Lᵣ:={x∈ℝⁿ: Ψ(x)≤r}. Choose any ε>0 such that ε∉Ψ(H⁻¹(0)). The weak Morse-Sard condition in Assumption 3.1 implies that ε can be chosen arbitrarily small, and Lemma B.9 implies that there are infinitely many k such that xₖ ∈L_ε. For any xₖ ∈L_ε, by the continuity of Ψ, we have that there exists α>0 such that dist(xₖ,ℝⁿ ∖L₂ε)>α. By (24), for sufficiently large k, we have that ‖ xₖ₊₁-xₖ ‖<α. Therefore, for all large k, xₖ ∈L_ε implies that xₖ₊₁ ∈L₂ε. Now, we define the last entrance and the first subsequent exit times,

$$k_{i}=\max\lbrace m\geq j_{i-1}:x_{m}\in\mathcal{L}_{\epsilon}\rbrace ,\quad j_{i}=\min\lbrace m\geq k_{i}:x_{m}\in\mathbb{R}^{n}\setminus\mathcal{L}_{2\epsilon}\rbrace .$$

We prove that such upcrossing occurs for finite times. Otherwise, if there exists {kᵢ} such that limᵢ→∞ kᵢ=∞, then Theorem 3.1 indicates that, up to a subsequence, ∇ϕ(x^λ_η(kᵢ)(·)) converges to ∇ϕ(x̅(·)), where x̅(·) is a trajectory of (6). By the definition of kᵢ, we have Ψ(x_kᵢ)≤ε, Ψ(x_kᵢ+1)>ε. By (24) and the continuity of $\Psi\circ\nabla\phi^{\ast }$, we have that limᵢ→∞ Ψ(x_kᵢ)=limᵢ→∞ Ψ(x_kᵢ+1)=ε. Recall that x^λ_η(kᵢ)(0)=x_kᵢ, therefore, Ψ(x̅(0))=limᵢ→∞ Ψ(x_kᵢ)=ε. Since x̅(0) is not in the stable set, there exists T>0, such that

$$\Psi(\bar{x}(T))\lt \sup_{s\in[0,T]}\Psi(\bar{x}(s))\leq\Psi(\bar{x}(0))=\epsilon.$$

Then, there exists δ>0, such that Ψ(x̅(T))≤ε-2δ. Moreover, for sufficiently large i, we have

$$\sup_{s\in[0,T]}\Psi(x^{\lambda_{\eta}(k_{i})}(s))\leq\sup_{s\in[0,T]}\Psi(\bar{x}(s))+\sup_{s\in[0,T]}|\Psi(x^{\lambda_{\eta}(k_{i})}(s))-\Psi(\bar{x}(s))|\leq 2\epsilon.$$

The last inequality comes from the uniform convergence of {∇ϕ(x^λ_η(kᵢ)(·))} in C(ℝ₊,ℝⁿ). This implies that for all large i, {x(λ_η(kᵢ)+s):s∈[0,T]} ⊂L₂ε. Thus it holds that λ_η(jᵢ)>λ_η(kᵢ)+T. Let lᵢ=max{m:λ_η(kᵢ)≤λ_η(m)≤λ_η(kᵢ)+T}. Then ‖ ∇ϕ(x_lᵢ)-∇ϕ(x^λ_η(kᵢ)(T))‖≤‖ ∇ϕ(x_lᵢ)-∇ϕ(x_lᵢ+1)‖ → 0, and hence ‖ ∇ϕ(x_lᵢ)-∇ϕ(x̅(T))‖ → 0 as i→∞. By the continuity of $\Psi\circ\nabla\phi^{\ast }$, we have that Ψ(x_lᵢ)≤ε-δ for all large i. By the definition of kᵢ and jᵢ, we have that λ_η(jᵢ)<λ_η(lᵢ)≤λ_η(kᵢ)+T, which leads to a contradiction. Therefore, for all large k, Ψ(xₖ)≤2ε. Since ε can be chosen arbitrarily small, it holds that limₖ→∞ Ψ(xₖ)=0. This completes the proof.

Appendix Appendix C: Proof of Theorem 3.3

Lemma C.11

Suppose Assumption 3.3 holds for the sequence of noises {ξₖ} and stepsizes {ηₖ}. Moreover, the Lyapunov function Ψ associated with H is coercive. Then for any ε>0 and T>0, there exists η_ε>0 such that for any {ηₖ} satisfying $\limsup_{k\to+\infty}\eta_{k}\leq\eta_{\varepsilon}$, almost surely, it holds that

$$\limsup_{s\to+\infty}\sup_{s\leq i\leq{\Lambda_{\eta}(\lambda_{\eta}(s)+T)}}\left\Vert \sum_{k=s}^{i}\eta_{k}\xi_{k}\right\Vert \leq\varepsilon.$$

(26)

Proof

From Assumption 3.3, it holds for all s≥0 and any i satisfying s≤i≤Λ_η(λ_η(s)+T) that

‖ Σₖ₌ₛⁱ ηₖ ξₖ ‖≤‖ Σₖ₌ₛ^N·⌈(s)/(N)⌉-1 ηₖ ξₖ ‖ +‖ Σₖ=N·⌈(s)/(N)⌉^N·⌊(i)/(N)⌋-1 ηₖ ξₖ ‖ +‖ Σₖ=N·⌊(i)/(N)⌋ⁱ ηₖ ξₖ ‖=‖ Σₖ₌ₛ^N·⌈(s)/(N)⌉-1 ηₖ ξₖ ‖ +‖ Σₖ=N·⌊(i)/(N)⌋ⁱ ηₖ ξₖ ‖ .

Let M_ξ be the uniform bound of {ξₖ}. Therefore, for any any ε>0, choosing η_ε=(ε)/(2NM_ξ) guarantees that

$$\limsup_{s\to+\infty}\sup_{s\leq i\leq\Lambda_{\eta}(\lambda_{\eta}(s)+T)}\left\Vert \sum_{k=s}^{i}\eta_{k}\xi_{k}\right\Vert \leq 2NM_{\xi}\limsup_{s\to+\infty,~{}s\leq i\leq\Lambda_{\eta}(\lambda_{\eta}(s)+T)}\eta_{i}\leq\varepsilon.$$

This completes the proof.

Proof of Theorem 3.3 . For the differential inclusion (6) and the update scheme (8), define zₖ=∇ϕ(xₖ) and z(t)=∇ϕ(x(t)). Then, we have

$$\dot{z}\in-\mathcal{H}(\nabla\phi^{\ast }(z))\text{and}z_{k+1}=z_{k}-\eta_{k}(d_{k}+\xi_{k}).$$

(28)

Note that $\mathcal{H}\circ\nabla\phi^{\ast }$ has a closed graph. By assumption, ϕ is μ-strongly convex and L-smooth, with some μ>0. Set δ̃ₖ:=max{δₖ,Lδₖ}. By L-smoothness of ϕ, We have

$$(\mathcal{H}\circ\nabla\phi^{\ast })^{\tilde{\delta}_{k}}(z_{k})=\mathcal{H}(\nabla\phi^{\ast }(z_{k}+\tilde{\delta}_{k}\mathbb{B}^{n}))+\tilde{\delta}_{k}\mathbb{B}^{n}\supset\mathcal{H}\left(\nabla\phi^{\ast }(z_{k})+\frac{\tilde{\delta}_{k}}{L}\mathbb{B}^{n})\right)+\tilde{\delta}_{k}\mathbb{B}^{n}\supset\mathcal{H}(\nabla\phi^{\ast }(z_{k})+\delta_{k}\mathbb{B}^{n})+\delta_{k}\mathbb{B}^{n}.$$

Note that dₖ ∈Hᵟₖ(xₖ)=H(∇ϕ∗(zₖ)+δₖ Bⁿ)+δₖ Bⁿ, it follows that $d_{k}\in(\mathcal{H}\circ\nabla\phi^{\ast })^{\tilde{\delta}_{k}}(z_{k})$. Lemma 4 and (xiao2023convergence, , Theorem 3.5) illustrate that for any ε̃>0, there exists ηₘₐₓ,T>0 such that whenever limsupₖ≥0 ηₖ≤ηₘₐₓ and {ξₖ} is (ε̃,T,{ηₖ})-controlled, then {zₖ} is uniformly bounded and

$$\limsup_{k\to+\infty}\;\mathrm{dist}\left(z_{k},\lbrace z\in\mathbb{R}^{n}:0\in\mathcal{H}(\nabla\phi^{\ast }(z))\rbrace \right)\leq{\tilde{\varepsilon}}.$$

Because ϕ is μ strongly convex, we have

$$\limsup_{k\to+\infty}\;\mathrm{dist}\left(x_{k},\lbrace x\in\mathbb{R}^{n}:0\in\mathcal{H}(x)\rbrace \right)\leq\limsup_{k\to+\infty}\;\frac{1}{\mu}\mathrm{dist}\left(z_{k},\lbrace x\in\mathbb{R}^{n}:0\in\mathcal{H}(\nabla\phi^{\ast }(z))\rbrace \right)\leq\frac{\tilde{\varepsilon}}{\mu}.$$

Since ε̃ is arbitrary, this completes the proof.

References

• (1)

  Ajanthan, T., Gupta, K., Torr, P., Hartley, R., Dokania, P.: Mirror descent view for neural network quantization. In: International Conference on Artificial Intelligence and Statistics, pp. 2809–2817. PMLR (2021)

• (2)

  Alvarez, F., Bolte, J., Brahic, O.: Hessian Riemannian gradient flows in convex programming. SIAM J. Control Optim. 43(2), 477–501 (2004)

• (3)

  Amid, E., Warmuth, M.K.: Reparameterizing mirror descent as gradient descent. Adv. Neural Inf. Process. Syst. 33, 8430–8439 (2020)

• (4)

  Azizan, N., Lale, S., Hassibi, B.: Stochastic mirror descent on overparameterized nonlinear models. IEEE Trans. Neural Netw. Learn. Syst. 33(12), 7717–7727 (2021)

• (5)

  Bauschke, H.H., Bolte, J., Teboulle, M.: A descent lemma beyond Lipschitz gradient continuity: first-order methods revisited and applications. Math. Oper. Res. 42(2), 330–348 (2017)

• (6)

  Beck, A., Teboulle, M.: Mirror descent and nonlinear projected subgradient methods for convex optimization. Oper. Res. Lett. 31(3), 167–175 (2003)

• (7)

  Benaïm, M., Hofbauer, J., Sorin, S.: Stochastic approximations and differential inclusions. SIAM J. Control Optim. 44(1), 328–348 (2005)

• (8)

  Benamou, J.D., Carlier, G., Cuturi, M., Nenna, L., Peyré, G.: Iterative Bregman projections for regularized transportation problems. SIAM J. Sci. Comput. 37(2), A1111–A1138 (2015)

• (9)

  Bernstein, J., Wang, Y.X., Azizzadenesheli, K., Anandkumar, A.: signsgd: Compressed optimisation for non-convex problems. In: International Conference on Machine Learning, pp. 560–569. PMLR (2018)

• (10)

  Bolte, J., Daniilidis, A., Lewis, A., Shiota, M.: Clarke subgradients of stratifiable functions. SIAM J. Optim. 18(2), 556–572 (2007)

• (11)

  Bolte, J., Le, T., Pauwels, E.: Subgradient sampling for nonsmooth nonconvex minimization. SIAM J. Optim. 33(4), 2542–2569 (2023)

• (12)

  Bolte, J., Pauwels, E.: Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning. Math. Program. 188, 19–51 (2021)

• (13)

  Bolte, J., Sabach, S., Teboulle, M., Vaisbourd, Y.: First order methods beyond convexity and lipschitz gradient continuity with applications to quadratic inverse problems. SIAM J. Optim. 28(3), 2131–2151 (2018)

• (14)

  Bolte, J., Teboulle, M.: Barrier operators and associated gradient-like dynamical systems for constrained minimization problems. SIAM J. Control Optim. 42(4), 1266–1292 (2003)

• (15)

  Borkar, V.S.: Stochastic approximation: a dynamical systems viewpoint, vol. 48. Springer (2009)

• (16)

  Bregman, L.M.: The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming. USSR Comput. Math. Math. Phys. 7(3), 200–217 (1967)

• (17)

  Castera, C., Bolte, J., Févotte, C., Pauwels, E.: An inertial newton algorithm for deep learning. J. Mach. Learn. Res. 22(1), 5977–6007 (2021)

• (18)

  Chu, H.T., Liang, L., Toh, K.C., Yang, L.: An efficient implementable inexact entropic proximal point algorithm for a class of linear programming problems. Comput. Optim. Appl. 85(1), 107–146 (2023)

• (19)

  Clarke, F.H.: Optimization and Nonsmooth Analysis. SIAM (1990)

• (20)

  Coste, M.: An introduction to o-minimal geometry. Istituti editoriali e poligrafici internazionali Pisa (2000)

• (21)

  Cutkosky, A., Mehta, H.: Momentum improves normalized sgd. In: International Conference on Machine Learning, pp. 2260–2268. PMLR (2020)

• (22)

  Davis, D., Drusvyatskiy, D., Kakade, S., Lee, J.D.: Stochastic subgradient method converges on tame functions. Found. Comput. Math. 20(1), 119–154 (2020)

• (23)

  Defazio, A.: Understanding the role of momentum in non-convex optimization: Practical insights from a Lyapunov analysis. arXiv preprint arXiv:2010.00406 (2020)

• (24)

  Ding, K., Li, J., Toh, K.C.: Nonconvex stochastic Bregman proximal gradient method with application to deep learning. J. Mach. Learn. Res. 26(39), 1–44 (2025)

• (25)

  Ding, K., Xiao, N., Toh, K.C.: Adam-family methods with decoupled weight decay in deep learning. Trans. Mach. Learn. Res. (2025). URL https://openreview.net/forum?id=xVEHiAZ7uR

• (26)

  Dragomir, R.A., Taylor, A.B., d’Aspremont, A., Bolte, J.: Optimal complexity and certification of Bregman first-order methods. Math. Program. pp. 1–43 (2021)

• (27)

  Van den Dries, L., Miller, C.: Geometric categories and o-minimal structures. Duke Math. J. (1996)

• (28)

  Duchi, J.C., Ruan, F.: Stochastic methods for composite and weakly convex optimization problems. SIAM J. Optim. 28(4), 3229–3259 (2018)

• (29)

  Ghai, U., Lu, Z., Hazan, E.: Non-convex online learning via algorithmic equivalence. Adv. Neural Inf. Process. Syst. 35, 22,161–22,172 (2022)

• (30)

  Gunasekar, S., Lee, J., Soudry, D., Srebro, N.: Characterizing implicit bias in terms of optimization geometry. In: International Conference on Machine Learning, pp. 1832–1841. PMLR (2018)

• (31)

  Guo, Z., Xu, Y., Yin, W., Jin, R., Yang, T.: A novel convergence analysis for algorithms of the adam family. arXiv preprint arXiv:2112.03459 (2021)

• (32)

  Gürbüzbalaban, M., Ruszczyński, A., Zhu, L.: A stochastic subgradient method for distributionally robust non-convex and non-smooth learning. J. Optim. Theory Appl. 194(3), 1014–1041 (2022)

• (33)

  Hanzely, F., Richtárik, P.: Fastest rates for stochastic mirror descent methods. Comput. Optim. Appl. 79, 717–766 (2021)

• (34)

  He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 770–778 (2016)

• (35)

  Hochreiter, S., Schmidhuber, J.: Long short-term memory. Neural Comput. 9(8), 1735–1780 (1997)

• (36)

  Jelassi, S., Li, Y.: Towards understanding how momentum improves generalization in deep learning, 2022. In: URL https://openreview.net/forum

• (37)

  Krizhevsky, A., Hinton, G., et al.: Learning multiple layers of features from tiny images (2009)

• (38)

  Lan, G.: Policy mirror descent for reinforcement learning: Linear convergence, new sampling complexity, and generalized problem classes. Math. Program. 198(1), 1059–1106 (2023)

• (39)

  Le, T.: Nonsmooth nonconvex stochastic heavy ball. J. Optim. Theory Appl. 201(2), 699–719 (2024)

• (40)

  Li, Y., Ju, C., Fang, E.X., Zhao, T.: Implicit regularization of Bregman proximal point algorithm and mirror descent on separable data. arXiv preprint arXiv:2108.06808 (2021)

• (41)

  Li, Z., Wang, T., Lee, J.D., Arora, S.: Implicit bias of gradient descent on reparametrized models: On equivalence to mirror descent. Adv. Neural Inf. Process. Syst. 35, 34,626–34,640 (2022)

• (42)

  Liu, Y., Gao, Y., Yin, W.: An improved analysis of stochastic gradient descent with momentum. Adv. Neural Inf. Process. Syst. 33, 18,261–18,271 (2020)

• (43)

  Lu, H., Freund, R.M., Nesterov, Y.: Relatively smooth convex optimization by first-order methods, and applications. SIAM J. Optim. 28(1), 333–354 (2018)

• (44)

  Marcus, M., Santorini, B., Marcinkiewicz, M.A.: Building a large annotated corpus of english: The penn treebank (1993)

• (45)

  Nemirovskij, A.S., Yudin, D.B.: Problem complexity and method efficiency in optimization. Wiley-Interscience (1983)

• (46)

  Polyak, B.T.: Some methods of speeding up the convergence of iteration methods. USSR Comput. Math. Math. Phys. 4(5), 1–17 (1964)

• (47)

  Ramezani-Kebrya, A., Antonakopoulos, K., Cevher, V., Khisti, A., Liang, B.: On the generalization of stochastic gradient descent with momentum. J. Mach. Learn. Res. 25(22), 1–56 (2024)

• (48)

  Rockafellar, R., Wets, R.J.B.: Variational Analysis. Springer Verlag, Heidelberg, Berlin, New York (1998)

• (49)

  Rockafellar, R.T.: Convex Analysis, vol. 11. Princeton University Press (1997)

• (50)

  Rudin, W.: Principles of mathematical analysis. McGraw-Hill, Inc., New York (1964)

• (51)

  Rudin, W.: Functional analysis 2nd ed. McGraw-Hill, Inc., New York (1991)

• (52)

  Ruszczyński, A.: Convergence of a stochastic subgradient method with averaging for nonsmooth nonconvex constrained optimization. Optim. Lett. 14(7), 1615–1625 (2020)

• (53)

  Ruszczynski, A.: A stochastic subgradient method for nonsmooth nonconvex multilevel composition optimization. SIAM J. Control Optim. 59(3), 2301–2320 (2021)

• (54)

  Sun, H., Gatmiry, K., Ahn, K., Azizan, N.: A unified approach to controlling implicit regularization via mirror descent. J. Mach. Learn. Res. 24(393), 1–58 (2023)

• (55)

  Takahashi, S., Takeda, A.: Approximate Bregman proximal gradient algorithm for relatively smooth nonconvex optimization. Comput. Optim. Appl. 90(1), 227–256 (2025)

• (56)

  Wu, F., Rebeschini, P.: A continuous-time mirror descent approach to sparse phase retrieval. Adv. Neural Inf. Process. Syst. 33, 20,192–20,203 (2020)

• (57)

  Xiao, N., Ding, K., Hu, X., Toh, K.C.: Developing Lagrangian-based Methods for Nonsmooth Nonconvex Optimization. arXiv preprint arXiv:2404.09438 (2024)

• (58)

  Xiao, N., Hu, X., Liu, X., Toh, K.C.: Adam-family methods for nonsmooth optimization with convergence guarantees. J. Mach. Learn. Res. 25(48), 1–53 (2024)

• (59)

  Xiao, N., Hu, X., Toh, K.C.: Convergence guarantees for stochastic subgradient methods in nonsmooth nonconvex optimization. arXiv preprint arXiv:2307.10053 (2023)

• (60)

  Yang, L., Toh, K.C.: Bregman proximal point algorithm revisited: A new inexact version and its inertial variant. SIAM J. Optim. 32(3), 1523–1554 (2022)

• (61)

  You, Y., Gitman, I., Ginsburg, B.: Large batch training of convolutional networks. arXiv preprint arXiv:1708.03888 (2017)

• (62)

  You, Y., Li, J., Reddi, S., Hseu, J., Kumar, S., Bhojanapalli, S., Song, X., Demmel, J., Keutzer, K., Hsieh, C.J.: Large batch optimization for deep learning: Training bert in 76 minutes. arXiv preprint arXiv:1904.00962 (2019)

• (63)

  Zhan, W., Cen, S., Huang, B., Chen, Y., Lee, J.D., Chi, Y.: Policy mirror descent for regularized reinforcement learning: A generalized framework with linear convergence. SIAM J. Optim. 33(2), 1061–1091 (2023)

• (64)

  Zhang, S., He, N.: On the convergence rate of stochastic mirror descent for nonsmooth nonconvex optimization. arXiv preprint arXiv:1806.04781 (2018)