Subspace correction methods, Block coordinate descent methods, Domain decomposition methods, Randomized methods, Convex optimization

1 Introduction↩︎

The main purpose of this paper is to develop an abstract framework for randomized subspace correction methods in convex optimization. While subspace correction methods [1], [2] generalize a broad class of iterative algorithms, convex optimization itself encompasses a wide range of applications, including nonlinear partial differential equations (PDEs), imaging, and data science. Consequently, the proposed framework is flexible and applicable to a very broad range of problems. Moreover, it accommodates highly general settings, such as weaker smoothness and convexity assumptions [3], [4].

Subspace correction methods [2] follow a divide-and-conquer strategy by decomposing the original problem into local subproblems defined on subspaces, which are solved independently. Many classical and modern iterative methods, including block relaxation, domain decomposition, and multigrid methods, can be viewed as instances of subspace correction methods. The theory has evolved over the past decades, covering both linear [2], [5] and nonlinear problems [1], [6], [7].

Block coordinate descent methods are prominent examples of subspace correction methods for convex optimization. These methods solve local subproblems restricted to blocks of coordinates, often via gradient or proximal steps. More general updates, like upper bound minimization [8], [9], are also possible. Their computational efficiency has led to widespread adoption. Key early results include [10], [11], and a comprehensive survey appears in [12]. Recent advances cover convergence of cyclic [13], [14], randomized [15], [16], accelerated [17], [18], and parallel variants [19]–[21]. Applications include deep neural network training [22], [23].

In numerical analysis, domain decomposition and multigrid methods are essential examples of subspace correction methods. Their convergence for smooth convex optimization was studied in [1], [24], and later extended to constrained and nonsmooth cases [7], [25]–[27]. These methods have been applied to various nonlinear variational problems, including PDEs [28]–[30], variational inequalities [26], [27], [31], elastoplasticity [6], and mathematical imaging [32]–[34].

Subspace correction methods are classified as parallel or successive depending on the order of subproblem updates [2], [35]. In parallel methods (additive Schwarz), all subproblems are solved concurrently; in successive methods (multiplicative Schwarz), subproblems are solved sequentially.

Randomizing the order of subproblem updates leads to randomized subspace correction methods [36], [37], which often exhibit better performance compared to fixed-order approaches. For quadratic optimization or linear problems, improved worst-case convergence rates under randomization are established in [37]. A notable result in [38] shows that the worst-case complexities of cyclic coordinate descent methods and randomized coordinate descent methods are \(\mathcal{O} (J^{4} \kappa \log(1/\epsilon) )\) and \(\mathcal{O} (J^{2} \kappa \log(1/\epsilon) )\), respectively, where \(J\) denotes the dimension of the variable, \(\epsilon\) is a prescribed tolerance for the energy error, and \(\kappa\) is a problem-dependent factor. More precisely, by constructing a specific quadratic objective, it was shown that both complexity bounds can be attained simultaneously. This demonstrates that cyclic coordinate descent methods can be up to \(\mathcal{O}(J^{2})\) times worse than randomized coordinate descent methods, thereby clarifying the role of randomization: when the optimal ordering of subspaces is unknown, randomization can yield better performance than certain fixed-order strategies. One may also refer to [39] for an analysis of how randomization averages antisymmetric terms in descent inequalities. This mechanism is beneficial for the development of accelerated methods [17], [18] and has motivated extensive research on randomized methods in convex optimization [16], [20], [40], [41]. In particular, [40] studies randomized subspace correction methods for convex optimization with linear local problems, which can be viewed as a preliminary work for the present paper.

This paper introduces an abstract framework for randomized subspace correction methods for convex composite optimization [41] on reflexive Banach spaces, accommodating diverse levels of smoothness and convexity [3], [4]. The framework unifies a wide range of decomposition strategies, including block partitioning [13], [15], [16] and overlapping domain decompositions commonly used in the numerical solution of PDEs [2], [29], [30], [35]. It supports both exact and inexact local solvers, encompassing methods such as coordinate descent, Bregman descent [42]–[44], and constraint decomposition [45], [46]. We establish convergence theorems that extend recent results, including those in [16], [36]. More precisely, by identifying a structural relationship between randomized and parallel subspace correction methods (see [37] for the linear case), we analyze the convergence of randomized subspace correction methods by leveraging existing analyses of parallel methods, such as [7], [47], [48]. Furthermore, we provide new analyses that yield sharper estimates under stronger, yet commonly adopted, assumptions; see, for example, [16].

In 1, we summarize the main results of this paper.

The remainder of this paper is organized as follows. In 2, we present an abstract framework of randomized subspace correction methods for convex optimization. In 3, we derive convergence theorems under various conditions on the target problem. In 4, we explain how the proposed framework relates to existing results. In 5, we summarize possible applications of the proposed framework from diverse fields of science and engineering. In 6, we provide numerical results for randomized subspace correction methods. Finally, in 7, we conclude the paper with some remarks.

2 Subspace correction methods↩︎

This section presents an abstract framework for randomized subspace correction methods for convex optimization. In particular, we show that the convergence analysis of randomized subspace correction methods for convex optimization can be carried out within the framework of parallel subspace correction methods [7], extending the analogy previously established for linear problems [36], [37]. The proposed framework is highly versatile, accommodating diverse space decomposition settings for the model problem, a broad range of smoothness and convexity levels in the objective functional, and various types of inexact local solvers.

Let \(V\) be a reflexive Banach space equipped with the norm \(\| \cdot \|\). Its topological dual is denoted by \(V^*\), and the duality pairing between \(V^*\) and \(V\) is written as \[\langle p, v \rangle = p(v), \quad p \in V^*,\;v \in V.\]

Throughout this paper, we adopt the convention \(0/0 = 0\) for arguments of \(\sup\) and \(0/0 = \infty\) for arguments of \(\inf\).

2.1 Space decomposition and subspace correction↩︎

We consider the following abstract convex optimization problem: \[\label{model} \min_{v \in V} \left\{ E(v) := F(v) + G(v) \right\},\tag{1}\] where \(F \colon V \to \mathbb{R}\) is a Gâteaux differentiable and convex functional, and \(G \colon V \to \overline{\mathbb{R}}\) is a proper, convex, and lower semicontinuous functional. The problem 1 is referred to as a composite optimization problem [41], as it involves a nonsmooth term \(G\) in addition to the smooth term \(F\). We further assume that the energy functional \(E\) is coercive, which guarantees the existence of a minimizer \(u \in V\) for the problem 1 .

We assume that the solution space \(V\) of 1 admits a space decomposition of the form \[\label{space95decomposition} V = \sum_{j=1}^J V_j,\tag{2}\] where each \(V_j\), \(j \in [J] = \{ 1, 2, \dots, J \}\), is a closed subspace of \(V\). The space decomposition 2 covers various algorithms, including block coordinate descent methods [15], [16], domain decomposition methods [7], [35], and multigrid methods [1], [49]. It is well known [5] that the space decomposition 2 satisfies the stable decomposition property. Namely, for any \(q \in [1, \infty)\), we have \[\label{norm95stable} \sup_{\| w \| = 1} \inf_{\sum_{j=1}^J w_j = w} \left( \sum_{j=1}^J \| w_j \|^q \right)^{\frac{1}{q}} < \infty,\tag{3}\] where \(w \in V\) and \(w_j \in V_j\).

Subspace correction methods involve solving local problems defined on subspaces \(\{ V_j \}_{j=1}^J\). For a given \(v \in V\), the optimal residual in a subspace \(V_j\) is obtained by solving the local minimization problem \[\label{local95exact} \min_{w_j \in V_j} E(v + w_j).\tag{4}\] Alternating minimization methods [13], [50] and certain domain decomposition methods (see, e.g., [1], [26], [29], [30]) fall into the category of subspace correction methods with exact local solvers as in 4 . In contrast, block coordinate descent methods typically solve the local problem 4 inexactly, often using a single iteration of gradient descent [13], [15], proximal descent [16], [18], or Bregman descent [42]–[44]. Some methods further employ surrogate techniques, where 4 is replaced by an approximate problem with lower computational complexity; see, e.g., [28], [45], [46].

To encompass all these methods, following [7], [48], we consider local problems of the form \[\label{local95inexact} \min_{w_j \in V_j} \left\{ E_j(w_j; v) := F_j(w_j; v) + G_j(w_j; v) \right\},\tag{5}\] where \(F_j(\cdot ; v) \colon V_j \to \mathbb{R}\) and \(G_j(\cdot ; v) \colon V_j \to \overline{\mathbb{R}}\) are convex functionals for each \(v \in V\). The functionals \(F_j(\cdot ; v)\) and \(G_j(\cdot ; v)\) serve as approximations to the exact local functionals \(F(v + \cdot)\) and \(G(v + \cdot)\) on \(V\), respectively. An example of 5 corresponding to a proximal descent step is presented in 1. Additional examples can be found in [7].

The abstract parallel subspace correction method for solving the convex optimization problem 1 , based on the space decomposition 2 and local solvers 5 , is presented in 1.

Another type of subspace correction method is the successive subspace correction method, in which the local problems in the subspaces are solved sequentially. In this paper, we focus on a particular variant known as the randomized subspace correction method, where the order of the local problems is chosen randomly; see 2.

2.2 Descent property↩︎

In what follows, we denote by \(d\) and \(d_j\) the Bregman divergences associated with \(F\) and \(F_j\), respectively: \[\begin{align} d(w; v) = F(v + w) - F(v) - \langle F'(v), w \rangle, \quad & v, w \in V, \\ d_j(w_j; v) = F_j(w_j; v) - F_j(0; v) - \langle F_j'(0; v), w_j \rangle, \quad & v \in V,\;w_j \in V_j. \end{align}\]

To ensure the convergence of the randomized subspace correction method, we adopt the assumptions on the local problem 5 summarized in 1. We note that 1 provides a more general framework than several recent works, as it extends the smooth settings in [47], [48] to the nonsmooth case, and employs a broader local stability assumption (see 1(c)) than the one used in [7].

The constant \(\rho\) defined in ?? is always greater than or equal to \(1\) as a consequence of 1(a,b). In the case of linear problems, one can verify that \(\rho = 2\) [48], which is consistent with [5], [35]. A nonlinear example where \(\rho > 1\) is provided in [47]. In 1, which is a nonsmooth extension of [48], we show that 1 ensures that solving each local problem leads to a decrease in the global energy.

As a corollary, the energy sequence generated by 2 decreases monotonically; see 1.

In 2, we present a refined version (cf. [48]) of the generalized additive Schwarz lemma for the composite optimization problem 1 , originally introduced in [7].

2 shows that, to analyze the convergence rate of the parallel subspace correction method (1), it suffices to estimate the following quantity [47], [48]: \[\label{Psi} \Psi(u^{(n)}) := \min_{w \in V} \left\{ \langle F'(u^{(n)}), w \rangle + \inf_{w = \sum_{j=1}^J w_j} \sum_{j=1}^J ( d_j + G_j )(w_j; u^{(n)}) \right\} - J G(u^{(n)}).\tag{11}\]

In the randomized subspace correction method (2), the update at each iteration is determined by a randomly chosen subspace. Consequently, the total energy after one step, \(E(u^{(n+1)})\), is a random variable depending on the sampling of the index \(j \in [J]\). To analyze its expected descent behavior, we consider the conditional expectation of the energy \(\mathbb{E} [ E (u^{(n+1)}) \mid u^{(n)} ]\), which represents the expected value of the energy at the next iteration given the current iterate \(u^{(n)}\). Here, the roman font \(E\) denotes the objective functional, while \(\mathbb{E}\) denotes expectation. In 1, we show that the conditional expectation \(\mathbb{E} [ E (u^{(n+1)}) \mid u^{(n)} ]\) can be estimated using 11 , indicating that its analysis can proceed along similar lines as that of the parallel method.

Such nonuniform sampling is particularly beneficial when the contributions of the \((d_j + G_j)\)-terms vary significantly across \(j \in [J]\). In this case, one may choose smaller \(\gamma_j\) for indices with larger contributions of \((d_j + G_j)\) to control \(\Psi_{\boldsymbol{\gamma}}(u^{(n)})\); see, e.g., [52] for specific sampling strategies and the resulting improvements.

3 Convergence theorems↩︎

In this section, we present convergence theorems for the randomized subspace correction method (2) under various conditions on the energy functional \(E\). We note that the convergence analysis of 2 given in this section relies on existing results for the parallel subspace correction method (1) [7], [47]. In particular, convergence rates of 1 have been established in [7] for general convex problems and in [7] for sharp (cf. 3) convex problems.

Given the initial iterate \(u^{(0)} \in \operatorname{dom} G\) of 2, we define \[\label{K950} K_0 = \{ v \in V : E (v) \leq E (u^{(0)}) \}, \quad R_0 = \sup_{v \in K_0} \| v - u \|.\tag{15}\] The convexity and coercivity of \(E\) imply that \(K_0\) is bounded and convex, and in particular, \(R_0 < \infty\). Moreover, by 1, the sequence \(\{ u^{(n)} \}\) generated by 2 remains entirely within \(K_0\).

1 implies that, to estimate the convergence rate of 2, it suffices to estimate \(\Psi(u^{(n)})\) defined in 11 . The form of \(\Psi(u^{(n)})\) naturally motivates the following stable decomposition assumption (cf. [7]).

Examples of stable decompositions satisfying 2 will be provided in [Sec:Related] [Sec:Applications]; see also [7]. A notable observation made in [48] is that, in the case of smooth problems, i.e., when \(G = 0\) in 1 and \(G_j = 0\), \(j \in [J]\), in 5 , 2 need not be assumed, but instead holds automatically under a mild smoothness condition on each \(d_j\); see [Prop:stable].

The following lemma provides a preliminary estimate for \(\Psi(u^{(n)})\) under 2. Although the proof follows a similar argument to that in [7], we include it here for completeness.

If one can derive suitable bounds for the term \[\label{Psi95tilde} \tilde{\Psi}(t; u^{(n)}) := t \langle F'(u^{(n)}), u - u^{(n)} \rangle + G\big((1-t)u^{(n)} + tu\big) - G(u^{(n)}), \quad t \in [0,1],\tag{17}\] then, in view of 3, one can proceed to obtain quantitative convergence bounds. Sharper estimates for \(\tilde{\Psi}(t; u^{(n)})\) can be established under stronger assumptions on \(F\) and \(G\). Accordingly, we consider several such cases in the following subsections.

3.1 General problems↩︎

Without imposing additional assumptions on \(F\) and \(G\), we can still obtain the following upper bound for 17 using the convexity of \(F\) and \(G\): \[\label{Psi95tilde95bound} \tilde{\Psi} (t; u^{(n)}) \leq - t ( E(u^{(n)}) - E(u)).\tag{19}\] By combining 1, 3, and 19 , we obtain the following convergence theorem for the randomized subspace correction method (2).

The generality of the assumptions in 2 enables a broad range of applications, particularly in scenarios where \(F\) exhibits a weaker level of smoothness than the standard smoothness condition [3]; see [47], [53] for concrete examples.

3.2 Sharp problems↩︎

Meanwhile, in many applications, the energy functional \(E\) satisfies the sharpness condition [4], summarized in 3, which is also known as the Hölderian error bound of the Łojasiewicz inequality [11], [54].

If we additionally assume that 3 holds, then we can derive the following improved convergence theorem for the randomized subspace correction method.

3.3 Strongly convex problems↩︎

Next, we consider a stronger condition than 3, stated in 4. Recall that a functional \(H \colon V \to \overline{\mathbb{R}}\) is said to be \(\mu\)-strongly convex on a convex set \(K \subset V\) if \[H ((1-t)v + tw) \leq (1-t) H(v) + t H(w) - t(1-t) \frac{\mu}{2} \| v - w \|^2, \quad t \in [0,1],\;v,w \in K.\] Note that 4 implies 3 with \(p = 2\).

Note that \(\mu_K\geq \mu_{F,K}\). Under 4, we have the following upper bound for 17 : \[\label{Psi95tilde95bound95strong} \begin{align} \tilde{\Psi} (t; u^{(n)}) &\leq E( (1-t) u^{(n)} + t u) - E(u^{(n)}) - \frac{\mu_{F, K_0}}{2} t^2 \| u - u^{(n)} \|^2 \\ &\leq -t ( E(u^{(n)}) - E(u)) + \left( \frac{\mu_{K950} - \mu_{F, K_0}}{2} t^2 - \frac{\mu_{K950}}{2} t \right) \| u - u^{(n)} \|^2, \end{align}\tag{27}\] where the first and second inequalities follow from the strong convexity of \(F\) and \(E\), respectively.

As 4 imposes a stronger condition than 3, the convergence rate established in 3(b) is guaranteed under 4. However, by using 27 , we can derive an even sharper estimate for the convergence rate, as presented in 4, when \(q = 2\) in 2. A similar result appears in [30].

4 Derivation of related methods↩︎

To highlight the versatility of the randomized subspace correction framework presented in this paper, we demonstrate in this section how several related results can be derived from it, including the randomized subspace correction method for linear problems [36], [37], block coordinate descent methods [15], [16], [18], and operator splitting methods [55].

Throughout this section, we assume that \(V\) is a Hilbert space equipped with an inner product \((\cdot, \cdot)\) and the induced norm \(\| \cdot \|\). Moreover, we identify \(V\) with its topological dual space \(V^*\) (cf. [49]).

4.1 Linear problems↩︎

The randomized subspace correction method for linear problems has been previously studied in [36], [37]. Here, we demonstrate how the general framework introduced in this paper recovers these existing results.

In 1 , we set \[F(v) = \frac{1}{2} (Av, v) - (f, v), \quad G(v) = 0, \quad v \in V,\] where \(A \colon V \to V\) is a linear operator induced by a continuous, symmetric, and coercive bilinear form on \(V\), and \(f \in V\). Then, it is readily observed that 1 reduces to the linear problem \[\label{linear} Au = f.\tag{31}\]

In the local problem 5 , for each \(j \in [J]\), we set \[\label{local95linear} F_j (w_j; v) = F(v) + (F'(v), w_j) + \frac{1}{2} (R_j^{-1} w_j, w_j), \quad G_j (w_j; v) = 0, \quad w_j \in V_j,\;v \in V,\tag{32}\] where \(R_j \colon V_j \to V_j\) is a linear operator induced by a continuous, symmetric, and coercive bilinear form on \(V_j\). Then, the solution \(\hat{w}_j\) of 5 is given by \[\hat{w}_j = - R_j Q_j (Av - f),\] where \(Q_j \colon V \to V_j\) denotes the orthogonal projection onto \(V_j\). Hence, we observe that 2 reduces to the randomized subspace correction method for linear problems introduced in [36], [37].

To analyze the algorithm, it suffices to verify 1 2 4. We adopt the \(\| \cdot \|_A\)-norm defined by \(\| \cdot \|_A := (A \cdot, \cdot)^{\frac{1}{2}}\). Note that 1(a, b) are trivially satisfied due to 32 . In 1(c), we have \(\rho = 2\) [48], and the assumption reduces to the following condition: for some \(\omega \in (0, 2)\), we require \[\label{local95linear95stable} (A w_j, w_j) \leq \omega (R_j^{-1} w_j, w_j), \quad w_j \in V_j,\tag{33}\] which corresponds to the standard assumption on local solvers, as found in, e.g., [35] and [5]. Moreover, [Prop:stable] ensures that 2 holds. Specifically, we have \(q = 2\), and for any bounded convex subset \(K\) of \(V\), it follows that \[C_K \leq \sup_{w \in V} \inf_{w = \sum_{j=1}^J w_j} \frac{\sum_{j=1}^J (R_j^{-1} w_j, w_j)}{(A w, w)} = \lambda_{\min}(T)^{-1},\] where the operator \(T \colon V \to V\) is given by \[\label{T} T = \sum_{j=1}^J R_j Q_j A,\tag{34}\] and the last equality follows from the well-known estimate in [5]; see also [7]. Finally, we observe that 4 holds with \[\mu_K = \mu_{F,K} \geq 1\] for any bounded convex subset \(K\) of \(V\). Therefore, under the condition 33 , the convergence estimate in [Rem:conv95strong] yields \[\mathbb{E}[\| u^{(n)} - u \|_A^2] \leq \left( 1 - \frac{\theta \lambda_{\min}(T)}{J} \right)^{n} \| u^{(0)} - u \|_A^2, \quad n \geq 0,\] where \(\theta\) and \(T\) were given in ?? and 34 , respectively (cf. [36]).

4.2 Block coordinate descent methods↩︎

Block coordinate descent methods are important instances of subspace correction methods for convex optimization. We discuss how randomized block coordinate descent methods [15], [16], [18] can be interpreted within the framework of the randomized subspace correction method; see also [55] and [7] for related discussions.

We begin by presenting the standard setting for block coordinate descent methods [16], [18]. Assume that the space \(V\) is given by the direct sum of subspaces \(V_j\), \(j \in [J]\): \[V = \bigoplus_{j=1}^J V_j,\] so that each \(v \in V\) can be represented in block form as [21] \[\label{block} v = (v_1, v_2, \dots, v_J), \quad v_j \in V_j.\tag{35}\]

In the composite optimization problem 1 , we assume that \(F\) is block smooth. That is, there exist positive constants \(L_j\), \(j \in [J]\), such that \[\label{F95RBCD} F(v + w_j) \leq F(v) + ( F'(v), w_j ) + \frac{L_j}{2} \| w_j \|^2, \quad v \in V,\;w_j \in V_j.\tag{36}\] In addition, we assume that the functional \(G\) is block separable, meaning it admits the decomposition \[\label{G95RBCD} G(v) = \sum_{j=1}^J G^j(v_j), \quad v \in V,\tag{37}\] for some functionals \(G^j \colon V_j \to \overline{\mathbb{R}}\).

In the local problem 5 , for each \(j \in [J]\), we set \[\label{local95RBCD} \begin{align} F_j(w_j; v) &= F(v) + ( F'(v), w_j ) + \frac{L_j}{2} \| w_j \|^2, \\ G_j(w_j; v) &= G ( v + w_j), \end{align} \quad w_j \in V_j,\;v \in V.\tag{38}\] Then, it is readily seen that the local problem 5 computes the increment of a single proximal descent step with step size \(1 / L_j\) applied to the \(j\)th block coordinate; see 1. Consequently, 2 reduces to the randomized block coordinate descent method considered in [16], [18].

Next, we demonstrate that our convergence theory recovers the results established in [16] for the randomized block coordinate descent method. It is straightforward to verify that 1 holds with \(\omega = 1\), based on 36 , 37 , and 38 . Moreover, using 37 , we obtain (cf. [7]) \[\sum_{j=1}^J G_j(w_j; v) = G(v + w) + (J - 1) G(v), \quad v,w \in V.\] From 36 and 38 , we have \[\sum_{j=1}^J d_j(w_j; v) \leq \frac{1}{2} \| w \|_L^2, \quad v, w \in V,\] where the norm \(\| \cdot \|_L\) is defined by \[\| v \|_L = \left( \sum_{j=1}^J L_j \| v_j \|^2 \right)^{\frac{1}{2}}, \quad v \in V.\] That is, 2 holds with the \(\| \cdot \|_L\)-norm, \(q = 2\), and \(C_K = 1\) for any bounded convex subset \(K \subset V\). Consequently, by 2, the randomized block coordinate descent method satisfies \[\mathbb{E}[E(u^{(n)})] - E(u) \leq \frac{2 J R_0^2}{n + 2 J R_0^2 / ( E(u^{(0)}) - E(u) )}, \quad n \geq 0,\] provided that \(E(u^{(0)}) - E(u)\) is sufficiently small. Here, \(R_0\) is defined by 15 with respect to the \(\| \cdot \|_L\)-norm. This result is consistent with [16] and [18].

Now suppose further that \(F\) and \(G\) are \(\mu_F\)- and \(\mu_G\)-strongly convex with respect to the \(\| \cdot \|_L\)-norm, respectively, for some \(\mu_F, \mu_G \geq 0\) with \(\mu_F + \mu_G > 0\). In this case, 4 implies the following linear convergence rate: \[\mathbb{E}[E(u^{(n)})] - E(u) \leq \left(1 - \frac{1}{J} \frac{\mu_F + \mu_G}{1 + \mu_G}\right)^{n} (E(u^{(0)}) - E(u)), \quad n \geq 0,\] which agrees with [16].

4.3 Operator splitting methods↩︎

By utilizing the duality between subspace correction and operator splitting methods developed in [55], we can derive a randomized operator splitting method from the randomized subspace correction method.

As a model problem for operator splitting methods, we consider the following optimization problem involving the sum of multiple convex functionals: \[\label{model95OS} \min_{v \in V} \left\{ F(v) + \sum_{j=1}^J G_j (B_j v) \right\},\tag{39}\] where \(V\) and each \(W_j\), \(j \in [J]\), are Hilbert spaces, \(B_j \colon V \to W_j\) are continuous linear operators, and \(F \colon V \to \mathbb{R}\) and \(G_j \colon W_j \to \overline{\mathbb{R}}\) are proper, convex, and lower semicontinuous functionals. For simplicity, we assume that \(F\) is strongly convex and smooth so that \(F'\) is invertible [56].

A randomized Peaceman–Rachford-type splitting method for solving 39 is presented in 3. At each iteration of 3, an index \(j \in [J]\) is selected at random, and an optimization problem involving only the functional \(G_j\) is solved. Recall that \(d\) denotes the Bregman divergence associated with \(F\).

Similar to [55], 5 shows that 3 is a dualization of the randomized subspace correction method applied to the following dual problem: \[\label{model95OS95dual} \min_{(p_j)_{j=1}^J \in \bigoplus_{j=1}^J W_j} \left\{ F^* \left( - \sum_{j=1}^J B_j^* p_j \right) + \sum_{j=1}^J G_j^* (p_j) \right\},\tag{40}\] where \(F^* \colon V \to \mathbb{R}\) and \(G_j^* \colon W_j \to \overline{\mathbb{R}}\) denote the Legendre–Fenchel conjugates of \(F\) and \(G_j\), respectively, and \(B_j^*\) denotes the adjoint of \(B_j\). Note that 40 is an instance of 1 in which the nonsmooth part is block separable (cf. 37 ).

5 Applications↩︎

In this section, we present applications of the randomized subspace correction method to a range of problems arising in diverse areas of science and engineering. See 2 for a concise summary.

Following standard convention in the literature (see, e.g., [2]), we write \(x \lesssim y\), or equivalently \(y \gtrsim x\), to mean that there exists a constant \(C>0\), independent of the relevant parameters, such that \(x \leq C y\). Furthermore, we write \(x \eqsim y\) when both \(x \lesssim y\) and \(x \gtrsim y\) hold.

5.1 Linear problems↩︎

A fundamental class of problems covered by the randomized subspace correction method is the class of linear problems [36], [37], as discussed in 4.1. Given the extensive literature on subspace correction methods for linear systems, particularly those arising from the numerical discretization of elliptic PDEs, we omit detailed discussion for brevity. We refer the reader to [35], [57], [58] and the references therein.

5.2 Nonlinear partial differential equations↩︎

The randomized subspace correction method is also applicable to nonlinear PDEs that admit convex variational formulations. As an example, we briefly present the \(s\)-Laplacian problem, which was also considered in [1], [29]. For other examples of nonlinear PDEs, one may refer to, e.g., [28], [30].

We consider the following nonlinear problem: \[\begin{align} - \operatorname{div} ( | \nabla u |^{s-2} \nabla u ) = f & \quad \text{ in } \Omega, \\ u = 0 & \quad \text{ on } \partial \Omega, \end{align}\] where \(\Omega \subset \mathbb{R}^2\) is a bounded polygonal domain, \(s \in (1, \infty)\) with \(s \neq 2\), and \(f \in W^{-1, s^*} (\Omega)\) with \(1/s + 1/s^* = 1\). It is well-known that the above problem admits a weak formulation given by the following convex optimization problem: \[\label{sLap} \min_{v \in W_0^{1,s} (\Omega) } \left\{ \frac{1}{s} \int_{\Omega} | \nabla v |^s \,dx - \langle f, v \rangle \right\}.\tag{41}\]

To numerically solve 41 , we employ a finite element discretization. Let \(\mathcal{T}_h\) be a quasi-uniform triangulation of \(\Omega\), where \(h\) denotes the characteristic element diameter. We denote by \(S_h(\Omega)\) the lowest-order Lagrangian finite element space defined on \(\mathcal{T}_h\), incorporating the homogeneous essential boundary condition. Then, the finite element approximation of 41 is given by \[\label{sLap95FEM} \min_{v \in S_h(\Omega)} \left\{ \frac{1}{s} \int_{\Omega} |\nabla v|^s \,dx - \langle f, v \rangle \right\}.\tag{42}\] We observe that 42 is an instance of the abstract problem 1 . Namely, 42 corresponds to 1 with \[V = S_h(\Omega), \quad F(v) = \frac{1}{s} \int_{\Omega} |\nabla v|^s \,dx - \langle f, v \rangle, \quad G(v) = 0.\]

For the space decomposition 2 , two-level overlapping Schwarz and multigrid decompositions with exact local solvers were studied in [7], [29] and [1], respectively. Here, we analyze the two-level overlapping Schwarz method for solving 42  [1], [7], [29] within the proposed framework. Let \(\mathcal{T}_H\) be a quasi-uniform triangulation of \(\Omega\) such that \(\mathcal{T}_h\) is a refinement of \(\mathcal{T}_H\), where \(H\) denotes the characteristic element diameter of \(\mathcal{T}_H\). Let \(\{ \Omega_i \}_{i=1}^N\) be a quasi-uniform overlapping domain decomposition of \(\Omega\), where each subdomain \(\Omega_i\) is a union of \(\mathcal{T}_h\)-elements with diameter of order \(H\). The overlap width is denoted by \(\delta\). We note that, for any bounded and convex subset \(K \subset V\), the following estimate holds [7]: \[\label{sLap95Bregman} \frac{\alpha_{K}}{p} \| u - v \|_{W^{1,s} (\Omega)}^p \leq F(u) - F(v) - \langle F'(v), u - v \rangle \leq \frac{\beta_{K}}{q} \|u - v \|_{W^{1,s} (\Omega)}^q, \quad u,v \in K,\tag{43}\] where \[p = \max \{ s, 2 \}, \quad q = \min \{ s, 2 \},\] and \(\alpha_K\), \(\beta_K\) are positive constants depending on \(K\).

In the two-level overlapping Schwarz method, we set \[\label{two95level} V_{N+1} = S_H (\Omega), \quad V_i = S_h (\Omega_i), \quad i \in [N],\tag{44}\] where \(S_H (\Omega)\) and \(S_h (\Omega_i)\) are defined analogously to \(S_h (\Omega)\). For the two-level space decomposition with \(J=N+1\), \[\label{two95level95space95decomposition} V = \sum_{j=1}^J V_j,\tag{45}\] the following stable decomposition result for 44 with respect to the \(\| \cdot \|_{W^{1,s}} (\Omega)\)-norm is available in [1], [7], [29]: for \(w \in V\), there exist \(w_j \in V_j\), \(1 \leq j \leq J\), such that \(w = \sum_{j=1}^J w_j\) and \[\label{sLap95stable} \sum_{j=1}^J \| w_j \|_{W^{1,s} (\Omega)}^s \lesssim \left( 1 + \left( \frac{H}{\delta} \right)^{s-1} \right) \| w \|_{W^{1,s} (\Omega)}^s.\tag{46}\] By 43 and 46 , 2 holds with \[C_{K950} \eqsim 1 + \left( \frac{H}{\delta} \right)^{\frac{q(s-1)}{s}}.\] Assuming that all local problems are solved exactly, i.e., using 4 , 1 holds trivially. Moreover, by 43 , 3 holds with \(\mu \eqsim 1\). Invoking 3, we obtain \[\mathbb{E}[E(u^{(n)})] - E(u) \lesssim \frac{1 + (H/\delta)^{\frac{q(s-1)}{s}}}{(n/J)^{\frac{p(q-1)}{p-q}}}\] That is, the two-level randomized Schwarz method for solving 42 converges sublinearly with rate \(\mathcal{O} (n^{- \frac{p(q-1)}{p-q}} )\).

5.3 Variational inequalities↩︎

Another important class of problems is variational inequalities, which find applications in computational mechanics and optimal control. As an illustrative example, we consider a second-order problem [27], [46]; see [6], [26], [31] for further examples.

We consider the following variational inequality: find \(u \in K\) such that \[\int_{\Omega} \nabla u \cdot \nabla v \,dx - \langle f, v \rangle \geq 0, \quad v \in K,\] where \(\Omega \subset \mathbb{R}^2\) is a bounded polygonal domain, \(f \in H^{-1}(\Omega)\), and \(K\) is a subset of \(H_0^1(\Omega)\) representing a pointwise inequality constraint: \[K = \{ v \in H_0^1(\Omega) : v \geq g \text{ a.e. in } \Omega \},\] for some \(g \in C(\Omega)\). This problem admits the equivalent optimization formulation \[\min_{v \in K} \left\{ \frac{1}{2} \int_{\Omega} |\nabla v|^2 \,dx - \langle f, v \rangle \right\}.\] Using a finite element discretization defined on \(S_h(\Omega)\), we obtain the discrete problem \[\label{VI95FEM} \min_{v \in S_h (\Omega)} \left\{ \frac{1}{2} \int_{\Omega} |\nabla v|^2 \,dx - \langle f, v \rangle + \chi_{K_h} (v) \right\},\tag{47}\] where \(K_h = K \cap S_h(\Omega)\), and \(\chi_{K_h}\) denotes the indicator functional of \(K_h\), defined as \(\chi_{K_h}(v) = 0\) if \(v \in K_h\) and \(\chi_{K_h}(v) = \infty\) otherwise. We observe that 47 is an instance of the general form 1 . Specifically, we obtain 47 if we set \[V = S_h(\Omega), \quad F(v) = \frac{1}{2} \int_{\Omega} |\nabla v|^2 \,dx - \langle f, v \rangle, \quad G(v) = \chi_{K_h}(v).\]

Two-level overlapping domain decomposition and multigrid decomposition with exact local solvers were studied in [25], [27]. In addition, constraint decomposition methods, based on localized constraints and interpretable as instances of inexact local solvers within our framework [7], were considered in [45], [46]. Here, we analyze the two-level overlapping Schwarz method for solving 47 , using the same subspaces as in 44 and the space decomposition 45 . We first observe that, since \(F\) is quadratic, by the Poincaré inequality we readily have \[\label{VI95Bregman} F(u) - F(v) - \langle F'(v), u-v \rangle \eqsim \| u - v \|_{H^1 (\Omega)}^2, \quad u,v \in V.\tag{48}\] Moreover, [30] (see also [27], [46]) proves the following stable decomposition result for 44 with respect to the \(\| \cdot \|_{H^1 (\Omega)}\)-norm: for \(v,w \in V\), there exists \(w_j \in V_j\), \(1 \leq j \leq J\), such that \(w = \sum_{j=1}^J w_j\) and \[\label{VI95stable} \begin{align} \sum_{j=1}^J \| w_j \|_{H^1 (\Omega)}^2 \lesssim \left(1 + \log \frac{H}{h} \right) \left( 1 + \frac{H}{\delta} \right) \| w \|_{H^1 (\Omega)}^2, \\ \sum_{j=1}^J G(v + w_j) \leq G(v + w) + J G(v). \end{align}\tag{49}\] By 48 and 49 , 2 holds with \(q=2\) and \[C_V \eqsim \left( 1 + \log \frac{H}{h} \right) \left(1 + \frac{H}{\delta} \right).\] Assuming that all local problems are solved exactly, i.e., using 4 , 1 holds trivially. Moreover, by 48 , 4 holds with \(\mu \eqsim 1\). Invoking 4, we obtain \[\mathbb{E} [E (u^{(n)})] - E(u) \leq \left( 1- \frac{1}{J C_V} \right)^n ( E (u^{(0)}) - E(u) ) .\] That is, the two-level randomized Schwarz method for solving 47 converges linearly.

5.4 Total variation minimization↩︎

Total variation minimization is a fundamental problem in mathematical imaging; see, e.g., [39]. Given a bounded polygonal domain \(\Omega \subset \mathbb{R}^2\), we consider the variational problem \[\label{TV} \min_{v \in BV(\Omega)} \left\{ \frac{1}{2} \int_{\Omega} (v - f)^2 \,dx + TV(v) \right\},\tag{50}\] where \(TV(v)\) denotes the total variation of \(v\), \(BV(\Omega)\) is the space of functions of bounded variation, and \(f \in L^2(\Omega)\).

Designing subspace correction methods for solving 50 is particularly challenging due to the nonseparable structure of the total variation term [61]. Indeed, it was shown in [60] that standard domain decomposition methods generally fail to satisfy the stable decomposition condition stated in 2.

One viable approach is to instead consider the dual formulation of 50  [62], which reads: \[\label{TV95dual} \min_{\mathbf{p} \in H_0 (\operatorname{div}; \Omega)} \frac{1}{2} \int_{\Omega} (\operatorname{div} \mathbf{p} + f)^2 \,dx \quad \text{subject to} \quad \| \mathbf{p} \|_{L^{\infty} (\Omega)} \leq 1.\tag{51}\] By employing the Raviart–Thomas finite element discretization introduced in [63], we obtain the discrete problem \[\label{TV95dual95FEM} \min_{\mathbf{p} \in \mathbf{S}_h (\Omega)} \left\{ \frac{1}{2} \int_{\Omega} (\operatorname{div} \mathbf{p} + f)^2 \,dx + \chi_{\mathbf{K}_h} (\mathbf{p}) \right\},\tag{52}\] where \(\mathbf{S}_h (\Omega)\) denotes the lowest-order Raviart–Thomas finite element space on the mesh \(\mathcal{T}_h\), and \(\mathbf{K}_h \subset \mathbf{S}_h (\Omega)\) is a convex set encoding the constraint \(\| \mathbf{p} \|_{L^{\infty} (\Omega)} \leq 1\). The discrete problem 52 fits the abstract formulation 1 , with the following identifications: \[V = \mathbf{S}_h (\Omega), \quad F(\mathbf{p}) = \frac{1}{2} \int_{\Omega} (\operatorname{div} \mathbf{p} + f)^2 \,dx, \quad G(\mathbf{p}) = \chi_{\mathbf{K}_h} (\mathbf{p}).\]

Schwarz-type domain decomposition methods for the dual problem 52 were analyzed in [61], and related constraint decomposition techniques were introduced in [32]; see also [33], [34] for domain decomposition methods for finite difference discretizations of 51 . Here, we analyze the one-level overlapping Schwarz method for solving 52 . In this method, we assume that the domain \(\Omega\) is decomposed into overlapping subdomains \(\{ \Omega_j \}_{j=1}^J\) as before. We set \[V_j = \mathbf{S}_h (\Omega_j), \quad j \in [J],\] where \(\mathbf{S}_h (\Omega_j)\) is defined analogously to \(\mathbf{S}_h (\Omega)\). Then we obtain the one-level space decomposition 2 .

We observe that \(F\) satisfies \[\label{TV95Bregman} F (\mathbf{p}) - F(\mathbf{q}) - \langle F'(\mathbf{q}), \mathbf{p} - \mathbf{q} \rangle = \frac{1}{2} \| \operatorname{div} ( \mathbf{p} - \mathbf{q} ) \|_{L^2 (\Omega)}^2, \quad \mathbf{p},\;\mathbf{q} \in V.\tag{53}\] It was shown in [61] that the following one-level stable decomposition with respect to the \(\| \cdot \|_{H (\operatorname{div}; \Omega)}\)-norm holds: for \(\mathbf{q}, \mathbf{r} \in V\), there exist \(\mathbf{r}_j \in V_j\), \(j \in [J]\), such that \(\mathbf{r} = \sum_{j=1}^J \mathbf{r}_j\) and \[\label{TV95stable} \begin{align} \sum_{j=1}^J \| \operatorname{div} \mathbf{r}_j \|_{L^2 (\Omega)}^2 \lesssim \frac{1}{\delta^2} \| \mathbf{r} \|_{H (\operatorname{div}; \Omega)}^2, \\ \sum_{j=1}^J G(\mathbf{q} + \mathbf{r}_j) \leq G(\mathbf{q} + \mathbf{r}) + (J-1) G(\mathbf{q}). \end{align}\tag{54}\] By 53 and 54 , 2 holds with \(q=2\) and \[C_V \eqsim \frac{1}{\delta^2}.\] Assuming that all local problems are solved exactly, 1 holds. Invoking 2, we obtain \[\mathbb{E}[E(\mathbf{p}^{(n)})] - E(\mathbf{p}) \lesssim \frac{\delta^{-2}}{n/J},\] where \(\mathbf{p} \in V\) denotes a solution of 52 . That is, the one-level randomized Schwarz method for solving 52 converges sublinearly with rate \(\mathcal{O}(n^{-1})\).

Meanwhile, to tackle the primal problem 50 directly, one may employ the operator splitting approach described in 4.3 to design domain decomposition methods. In this context, it was shown in [55] that the resulting methods coincide with those proposed in [60]; see also [59].

5.5 Multinomial logistic regression↩︎

Multinomial logistic regression is a fundamental tool for classification problems in machine learning. Given a labeled dataset \(\{ (x_{i}, y_{i} )\}_{i=1}^{N} \subset \mathbb{R}^d \times [k]\), where each \(x_{i} \in \mathbb{R}^d\) represents a data point and \(y_{i} \in [k]\) its corresponding class label, the regression model is formulated as follows: \[\label{logistic95regression} \min_{\theta \in \mathbb{R}^{(d+1)k}} \left\{ \frac{1}{N} \left( \sum_{i=1}^{N} \operatorname{LSE}_k (X_{i}^T \theta) - \hat{x}^T \theta \right) + \frac{\alpha}{2} \| \theta \|^2 \right\},\tag{55}\] where \(\theta = [w_1^T, b_1, \dots, w_k^T, b_k]^T \in \mathbb{R}^{(d+1)k}\) is the parameter vector, \(\operatorname{LSE}_k\) denotes the log-sum-exp function over \(k\) classes, \(X_{i} \in \mathbb{R}^{(d+1)k \times k}\) is defined as \(X_{i} = I_k \otimes [x_{i}^T, 1]^T\), and \(\alpha\) is a positive regularization hyperparameter. The vector \(\hat{x} \in \mathbb{R}^{(d+1)k}\) is a constant vector depending only on the dataset; see [55] for details. Note that 55 is an instance of 39 with \(J = N\).

In big data settings, where the number of data points \(N\) is very large, operator splitting methods such as stochastic gradient descent [64], which process individual data points or mini-batches, are widely used to solve 55 . The randomized Peaceman–Rachford splitting algorithm presented in 3 is one such method.

Thanks to 5, the convergence analysis of 3 for solving 55 reduces to analyzing the randomized subspace correction method for the following dual problem [65], [66]: \[\label{logistic95regression95dual} \min_{(p_{i})_{i=1}^{N} \in (\mathbb{R}^{k})^N} \left\{ \frac{1}{2 N \alpha} \left\| \sum_{i=1}^{N} X_{i} p_{i} - \hat{x} \right\|^2 + \sum_{i=1}^{N} \operatorname{LSE}_k^* (p_{i}) \right\}.\tag{56}\] More precisely, we consider the randomized subspace correction method for solving 56 with a nonoverlapping \(N\)-block decomposition, i.e., in 2 we set \[J = N, \quad V = (\mathbb{R}^{k})^N, \quad V_i = \mathbb{R}^k, \quad i \in [N].\] Since the decomposition is nonoverlapping, each \(p \in V\) admits a unique decomposition \(p = (p_i)_{i=1}^N\) with respect to the direct sum \(\bigoplus_{i=1}^N V_i\). Moreover, 56 is an instance of 1 with \[F(p) = \frac{1}{2 N \alpha} \left\| \sum_{i=1}^{N} X_{i} p_{i} - \hat{x} \right\|^2, \quad G(p) = \sum_{i=1}^{N} \operatorname{LSE}_k^* (p_{i}).\] Since the term \(\sum_{i=1}^{N} \operatorname{LSE}_k^* (p_{i})\) is block separable, an argument analogous to that in 4.2 shows that 2 holds with \(q = 2\). Consequently, 2 applies to establish convergence of the randomized subspace correction method for solving 56 .

6 Numerical results↩︎

In this section, we present numerical results for some applications of the randomized subspace correction method. All algorithms were implemented in MATLAB R2024b and run on a desktop equipped with an AMD Ryzen 5 5600X CPU (3.7 GHz, 6C), 40 GB RAM, and Windows 10 Pro.

6.1 \(s\)-Laplacian problem↩︎

As the first example, we consider two-level overlapping Schwarz methods for the \(s\)-Laplacian problem 41 . The additive and randomized Schwarz methods correspond to the parallel and randomized subspace correction methods in 1 2, respectively.

In our numerical experiments, we adopt the setting in [29]. We set \(\Omega = (0,1)^2 \subset \mathbb{R}^2\), \(s \in \{1.5, 5\}\), and \(f = 1\). For the finite element discretization 42 and the two-level space decomposition 45 , we partition \(\Omega\) into \(2 \times 1/H \times 1/H\) uniform triangles to obtain a coarse triangulation \(\mathcal{T}_H\), and then refine it to obtain a fine triangulation \(\mathcal{T}_h\) with \(2 \times 1/h \times 1/h\) triangles. Each nonoverlapping subdomain \(D_i\), \(i \in [N]\) with \(N = 1/H \times 1/H\), is defined as a rectangular region consisting of two coarse triangles sharing a diagonal edge. The corresponding overlapping subdomain \(D_i'\) is obtained by enlarging \(D_i\) with layers of fine triangles of width \(\delta\). The decomposition \(\{ D_i' \}_{i=1}^N\) can be colored with four colors when \(\delta\) is sufficiently small; see [1]. Let \(\Omega_j\), \(j \in [N_c]\) with \(N_c = 4\), denote the union of all \(D_i'\) with the \(j\)th color. Defining \(V_j\), \(j \in [J]\) with \(J = N_c + 1 = 5\), as in 44 , we obtain the two-level space decomposition 45 . This coloring-based space decomposition enables parallel computation in 2: all subdomain problems in the same color can be solved in parallel, while in 1 all subdomain problems regardless of color can be solved in parallel. In 1, we choose the step size \(\tau = 1/5\) (cf. [29]).

In both 1 2, we initialize \(u^{(0)}\) randomly. For solving the local and coarse problems on \(V_j\), \(1 \le j \le J\), we use the adaptive Newton method in [67] with the stopping criterion \[\label{stop} \left| \frac{E_j^{(n)}(w_j^{(m+1)}) - E_j^{(n)}(w_j^{(m)})}{E_j^{(n)}(w_j^{(m+1)})} \right| < 10^{-12}.\tag{57}\] Here \(E_j^{(n)}\) denotes the energy functional associated with \(V_j\) at the \(n\)th outer iteration, and \(m\) is the number of inner iterations. To obtain a reference solution \(u_h\), we apply the adaptive Newton method [67] to the full problem 42 with sufficiently many iterations.

In 4, we depict convergence curves of the relative energy error \[\label{rel95energy95error} \frac{E(u^{(n)}) - E(u)}{E(u^{(0)}) - E(u)}\tag{58}\] with respect to the number of color solves generated by 1 2 based on the two-level space decomposition 45 . Note that each outer iteration of 1 is counted as \(J+1 = 5\) color solves. We observe that 2 converges faster than 1, reflecting the well-known fact that successive methods such as Gauss–Seidel converge faster than their parallel counterparts such as Jacobi. Moreover, 2 is robust with respect to the sampling of subspace indices: even the worst case over all experiments, represented by the upper boundary of the shaded region in 4, outperforms 1.

6.2 \(L^1\)-penalized problem↩︎

As the second example, we consider two-level overlapping Schwarz methods [30] for the \(L^1\)-penalized problem introduced in [68]: \[\label{L1} \min_{v \in H_0^1 (\Omega)} \int_{\Omega} \left( \frac{1}{2} | \nabla v |^2 - gv + \alpha | v | \right) \,dx,\tag{59}\] where \(\alpha > 0\). In our numerical experiments, we adopt the setting in [30]. We set \(\Omega = (0,1)^2 \subset \mathbb{R}^2\), \(g(x,y) = 10^3 x(1-x) \sin(\pi y)\), and \(\alpha \in \{10, 30\}\). For discretization and domain decomposition, we use the same settings as for the \(s\)-Laplacian problem 41 . Note that the corresponding finite element discretization of 59 is an instance of 1 with \[V = S_h (\Omega), \quad F(v) = \frac{1}{2} \int_{\Omega} | \nabla v |^2 \,dx - \int_{\Omega} g v \,dx, \quad G(v) = \alpha \int_{\Omega} | v | \,dx.\]

Since 59 is nonsmooth, the adaptive Newton method introduced above is not applicable. We employ the fast gradient method with adaptive restart (AFGM) [69]. For solving the local problems, we use AFGM with the stopping criterion 57 . For the coarse problems, we apply AFGM to the dual formulations of the coarse problems, akin to the approach introduced in [31]. A reference solution \(u_h\) is obtained by running AFGM on the full problem with sufficiently many iterations.

In 5, we depict convergence curves of 1 2 based on the two-level space decomposition 45 . A similar discussion applies as in the case of the \(s\)-Laplacian problem. Namely, 2 converges faster than 1 with respect to the number of color solves, and is also robust with respect to the sampling of subspace indices.

7 Conclusion↩︎

In this paper, we introduced an abstract framework for randomized subspace correction methods for convex optimization. We established convergence theorems that are applicable to a broad range of scenarios involving space decomposition, local solvers, and varying levels of problem regularity. Furthermore, we demonstrated that these theorems unify and extend several relevant recent results.

This work opens several avenues for future research. One important direction is the development of convergence theory for cyclic successive subspace correction methods. This is particularly relevant for multigrid methods [1], [25], where the hierarchical structure of subspaces plays a critical role. We note that a sharp convergence analysis of cyclic methods for linear problems can be found in [5], [70]. The block coordinate descent method with \(q = 2\) was analyzed in [13]; however, the results there are not sharp in the sense that the analysis does not recover the sharp estimate established for the linear case.

Another promising direction is to extend the proposed framework to nonconvex problems. Recent work [71] has shown that randomized block coordinate descent methods are effective for a certain class of nonconvex problems. Generalizing the randomized subspace correction framework to accommodate such settings remains an interesting and challenging problem.

References↩︎

---

1. Applied Mathematics and Computational Sciences Program, Computer, Electrical and Mathematical Science and Engineering Division, King Abdullah University of Science and Technology (KAUST), Thuwal 23955, Saudi Arabia (, , ).↩︎

2. Submitted to arXiv.↩︎