Parallel Douglas–Rachford Algorithms

• Parallel Douglas–Rachford algorithms are splitting methods that extend the classical two-operator framework to multi-operator settings, enabling blockwise resolvent evaluations.
• They reformulate complex problems into a product space with a diagonal constraint, effectively decoupling global tasks into parallel local updates.
• Applications include convex minimization, feasibility, and primal–dual formulations, with convergence guarantees based on monotone operator theory and fixed-point analysis.

Parallel Douglas–Rachford type algorithms are splitting methods that extend the classical two-operator Douglas–Rachford paradigm to sums of many operators, many functions, or coupled primal–dual systems. Their characteristic device is a reformulation in a product space together with a diagonal constraint, a normal-cone term, or a primal–dual lifting, so that the difficult global problem is replaced by blockwise resolvent evaluations and a consensus-enforcing averaging or projection step. Within this framework, the literature covers monotone inclusions with composite and parallel-sum structure, convex minimization of finitely many terms, convex feasibility, inconsistent models with normal solutions, and geodesically convex optimization on Hadamard manifolds (Bot et al., 2012, Alcantara et al., 6 Jan 2025, Bauschke et al., 2019, Censor et al., 2015, Bergmann et al., 2015, Sahu et al., 28 Sep 2025).

1. Product-space reformulation and the diagonal principle

A standard formulation starts from an $m$-operator inclusion

$0 \in (A_1+\cdots+A_m)(x)$

on a real Hilbert space $H$, and lifts it to $\mathcal H=H^m$. One defines the block-diagonal operator

$F(X)=A_1(x_1)\times\cdots\times A_m(x_m),$

the diagonal subspace

$D=\{X\in H^m:x_1=\cdots=x_m\},$

and the normal-cone operator $G=N_D$. Then

$$0\in(F+G)(X) \quad\Longleftrightarrow\quad x_1=\cdots=x_m=x \text{ and } 0\in\sum_i A_i(x).$$

The projector onto $D$ is explicit: $$P_D(X)=(\bar y,\bar y,\dots,\bar y), \qquad \bar y=\frac1m\sum_{i=1}^m x_i.$$ This is the basic Pierra/Campoy-style product-space reformulation used for multioperator Douglas–Rachford splitting (Alcantara et al., 6 Jan 2025).

For convex minimization of finitely many functions,

$0 \in (A_1+\cdots+A_m)(x)$0

the same principle appears in the product space $0 \in (A_1+\cdots+A_m)(x)$1, with

$0 \in (A_1+\cdots+A_m)(x)$2

so that the problem becomes

$0 \in (A_1+\cdots+A_m)(x)$3

The diagonal condition is therefore encoded as a single linear subspace constraint, and the prox of $0 \in (A_1+\cdots+A_m)(x)$4 splits into the $0 \in (A_1+\cdots+A_m)(x)$5 independent prox operators $0 \in (A_1+\cdots+A_m)(x)$6 (Bauschke et al., 2019).

On manifolds the same architecture persists. For $0 \in (A_1+\cdots+A_m)(x)$7 on a symmetric Hadamard manifold, one lifts to $0 \in (A_1+\cdots+A_m)(x)$8 and adds the diagonal set $0 \in (A_1+\cdots+A_m)(x)$9. In Hadamard-manifold formulations, the diagonal projection is expressed through a Riemannian centroid,

$H$0

and $H$1 (Bergmann et al., 2015, Sahu et al., 28 Sep 2025).

The significance of these constructions is structural rather than cosmetic. They convert a many-term problem into a two-term Douglas–Rachford problem in which one resolvent is block-separable and the other is an averaging or projection operator. This is the source of the parallelism that characterizes the entire class (Alcantara et al., 6 Jan 2025, Bauschke et al., 2019).

2. Canonical block-parallel Douglas–Rachford iterations in Hilbert spaces

For maximal $H$2-monotone operators $H$3, the block-parallel Douglas–Rachford algorithm on $H$4 takes the standard two-operator form

$H$5

Because $H$6 and $H$7 are separable, the iteration can be written as three explicit steps: $H$8

$H$9

$\mathcal H=H^m$0

The “shadow” sequence $\mathcal H=H^m$1 and the projected sequence $\mathcal H=H^m$2 both serve as candidate solution sequences (Alcantara et al., 6 Jan 2025).

The convergence theory is divided into regimes by the aggregate monotonicity parameter $\mathcal H=H^m$3. In the convex/monotone case, where all $\mathcal H=H^m$4, any $\mathcal H=H^m$5 and any $\mathcal H=H^m$6 yield the classical Douglas–Rachford behavior: $\mathcal H=H^m$7, and $\mathcal H=H^m$8 with $\mathcal H=H^m$9. If $F(X)=A_1(x_1)\times\cdots\times A_m(x_m),$0, then one chooses $F(X)=A_1(x_1)\times\cdots\times A_m(x_m),$1 small enough so that $F(X)=A_1(x_1)\times\cdots\times A_m(x_m),$2; in that case the method converges strongly, $F(X)=A_1(x_1)\times\cdots\times A_m(x_m),$3 is unique, and both $F(X)=A_1(x_1)\times\cdots\times A_m(x_m),$4 and $F(X)=A_1(x_1)\times\cdots\times A_m(x_m),$5 converge to it. In the finite-dimensional nonconvex setting considered in the same work, where $F(X)=A_1(x_1)\times\cdots\times A_m(x_m),$6 operators are Lipschitz gradients of $F(X)=A_1(x_1)\times\cdots\times A_m(x_m),$7 functions and the remaining one is a proper closed subdifferential, blockwise choices $F(X)=A_1(x_1)\times\cdots\times A_m(x_m),$8 and $F(X)=A_1(x_1)\times\cdots\times A_m(x_m),$9 guarantee that every cluster point of the shadow sequence is a stationary point, with residual $D=\{X\in H^m:x_1=\cdots=x_m\},$0 (Alcantara et al., 6 Jan 2025).

The implementation pattern is explicit: each sweep consists of $D=\{X\in H^m:x_1=\cdots=x_m\},$1 parallel resolvent calls, one global sum/average, and $D=\{X\in H^m:x_1=\cdots=x_m\},$2 parallel local updates, with no further synchronization needed until the next iteration. In this form, “parallel Douglas–Rachford” is not merely a variant of the two-set algorithm but a general multioperator template built around blockwise separability and a single consensus operation (Alcantara et al., 6 Jan 2025).

3. Primal–dual Douglas–Rachford type splittings for composite and parallel-sum inclusions

A distinct but closely related line studies primal–dual splittings for coupled monotone inclusions. In real Hilbert spaces $D=\{X\in H^m:x_1=\cdots=x_m\},$3 and $D=\{X\in H^m:x_1=\cdots=x_m\},$4, with maximally monotone operators $D=\{X\in H^m:x_1=\cdots=x_m\},$5, $D=\{X\in H^m:x_1=\cdots=x_m\},$6, bounded linear operators $D=\{X\in H^m:x_1=\cdots=x_m\},$7, and data $D=\{X\in H^m:x_1=\cdots=x_m\},$8, $D=\{X\in H^m:x_1=\cdots=x_m\},$9, the primal inclusion is formulated as

$G=N_D$0

and equivalently, writing $G=N_D$1 for the parallel sum of $G=N_D$2 and $G=N_D$3,

$G=N_D$4

The associated dual inclusion is

$G=N_D$5

Under the qualification

$G=N_D$6

the set of primal–dual solutions is nonempty (Bot et al., 2012).

The paper develops two inexact Douglas–Rachford type primal–dual algorithms. Algorithm 1 is a two-step Douglas–Rachford splitting obtained by lifting $G=N_D$7 to a product space $G=N_D$8 and applying an inexact DR iteration in a metric induced by a strongly positive self-adjoint operator $G=N_D$9. With step sizes $$0\in(F+G)(X) \quad\Longleftrightarrow\quad x_1=\cdots=x_m=x \text{ and } 0\in\sum_i A_i(x).$$0 and $$0\in(F+G)(X) \quad\Longleftrightarrow\quad x_1=\cdots=x_m=x \text{ and } 0\in\sum_i A_i(x).$$1 satisfying

$$0\in(F+G)(X) \quad\Longleftrightarrow\quad x_1=\cdots=x_m=x \text{ and } 0\in\sum_i A_i(x).$$2

relaxation parameters $$0\in(F+G)(X) \quad\Longleftrightarrow\quad x_1=\cdots=x_m=x \text{ and } 0\in\sum_i A_i(x).$$3, and summable error sequences, the method uses the resolvents $$0\in(F+G)(X) \quad\Longleftrightarrow\quad x_1=\cdots=x_m=x \text{ and } 0\in\sum_i A_i(x).$$4, $$0\in(F+G)(X) \quad\Longleftrightarrow\quad x_1=\cdots=x_m=x \text{ and } 0\in\sum_i A_i(x).$$5, and $$0\in(F+G)(X) \quad\Longleftrightarrow\quad x_1=\cdots=x_m=x \text{ and } 0\in\sum_i A_i(x).$$6; all computations for $$0\in(F+G)(X) \quad\Longleftrightarrow\quad x_1=\cdots=x_m=x \text{ and } 0\in\sum_i A_i(x).$$7 can be carried out in parallel. Algorithm 2 is a single-step Douglas–Rachford splitting in which each evaluation of $$0\in(F+G)(X) \quad\Longleftrightarrow\quad x_1=\cdots=x_m=x \text{ and } 0\in\sum_i A_i(x).$$8 and $$0\in(F+G)(X) \quad\Longleftrightarrow\quad x_1=\cdots=x_m=x \text{ and } 0\in\sum_i A_i(x).$$9 is reduced to a single forward–backward pass. Its price is a stricter condition,

$D$0

together with an auxiliary variable $D$1 (Bot et al., 2012).

The convergence results mirror the operator-theoretic structure. For Algorithm 1, if $D$2 and $D$3 are merely maximally monotone, then $D$4 converges weakly to a primal–dual solution. In finite dimensions, if $D$5 is bounded away from $D$6 and $D$7 are uniformly monotone, then $D$8 and $D$9 strongly. For Algorithm 2, weak convergence holds under the same range condition, while strong convergence in finite dimensions requires uniform monotonicity of $$P_D(X)=(\bar y,\bar y,\dots,\bar y), \qquad \bar y=\frac1m\sum_{i=1}^m x_i.$$0, $$P_D(X)=(\bar y,\bar y,\dots,\bar y), \qquad \bar y=\frac1m\sum_{i=1}^m x_i.$$1, and $$P_D(X)=(\bar y,\bar y,\dots,\bar y), \qquad \bar y=\frac1m\sum_{i=1}^m x_i.$$2 (Bot et al., 2012).

The practical comparison between the two schemes is explicit. Algorithm 1 uses two forward evaluations of $$P_D(X)=(\bar y,\bar y,\dots,\bar y), \qquad \bar y=\frac1m\sum_{i=1}^m x_i.$$3 and $$P_D(X)=(\bar y,\bar y,\dots,\bar y), \qquad \bar y=\frac1m\sum_{i=1}^m x_i.$$4 per iteration but allows the larger stepsize bound $$P_D(X)=(\bar y,\bar y,\dots,\bar y), \qquad \bar y=\frac1m\sum_{i=1}^m x_i.$$5. Algorithm 2 halves the matrix-vector products at the cost of the stricter bound $$P_D(X)=(\bar y,\bar y,\dots,\bar y), \qquad \bar y=\frac1m\sum_{i=1}^m x_i.$$6 and the additional auxiliary variable $$P_D(X)=(\bar y,\bar y,\dots,\bar y), \qquad \bar y=\frac1m\sum_{i=1}^m x_i.$$7. Both treat the resolvents in parallel and admit inexact computations via summable errors (Bot et al., 2012).

4. Multiple-set, convex-feasibility, and inconsistent convex formulations

Parallel Douglas–Rachford type algorithms also arise in convex feasibility and in convex minimization over finitely many terms. For

$$P_D(X)=(\bar y,\bar y,\dots,\bar y), \qquad \bar y=\frac1m\sum_{i=1}^m x_i.$$8

the product-space reformulation with diagonal subspace $$P_D(X)=(\bar y,\bar y,\dots,\bar y), \qquad \bar y=\frac1m\sum_{i=1}^m x_i.$$9 yields the parallel componentwise update

$0 \in (A_1+\cdots+A_m)(x)$00

This is precisely the Douglas–Rachford operator $0 \in (A_1+\cdots+A_m)(x)$01 written in coordinates, so each iteration consists of $0 \in (A_1+\cdots+A_m)(x)$02 proximal steps in parallel plus one averaging step. In the finite-dimensional theory developed for this model, under proper convex lsc assumptions, a constraint qualification $0 \in (A_1+\cdots+A_m)(x)$03, and nonemptiness of the normal solution set, the shadow sequence $0 \in (A_1+\cdots+A_m)(x)$04 converges to a limit $0 \in (A_1+\cdots+A_m)(x)$05 solving

$0 \in (A_1+\cdots+A_m)(x)$06

where $0 \in (A_1+\cdots+A_m)(x)$07. This is the mechanism that allows the method to remain meaningful even in possibly inconsistent cases (Bauschke et al., 2019).

For convex feasibility with closed convex sets $0 \in (A_1+\cdots+A_m)(x)$08, new algorithmic structures embed the basic two-set Douglas–Rachford operator $0 \in (A_1+\cdots+A_m)(x)$09 into larger parallel architectures. In String-Averaging DR, one partitions the index set into strings $0 \in (A_1+\cdots+A_m)(x)$10, composes two-set DR operators along each string, and averages the resulting string operators with weights $0 \in (A_1+\cdots+A_m)(x)$11. If $0 \in (A_1+\cdots+A_m)(x)$12, the sequence converges strongly to a point in $0 \in (A_1+\cdots+A_m)(x)$13. In Block-Iterative DR, one block is selected per iteration, parallel pairwise DR operators are computed inside that block, and their weighted average defines the next iterate; with $0 \in (A_1+\cdots+A_m)(x)$14, the swept iterates converge weakly, and strong convergence follows under the interior-point condition. The same work also defines an $0 \in (A_1+\cdots+A_m)(x)$15-set DR operator

$0 \in (A_1+\cdots+A_m)(x)$16

where $0 \in (A_1+\cdots+A_m)(x)$17 is a composite reflection, and studies mixtures

$0 \in (A_1+\cdots+A_m)(x)$18

These schemes include the cyclic Douglas–Rachford algorithm, the simultaneous DR operator, and averaged DR as special cases (Censor et al., 2015).

The operator-theoretic language in this literature is strongly quasi-nonexpansive and firmly nonexpansive. Compositions and convex combinations preserve these classes, which makes fixed-point analysis natural. The practical message is that “parallel DR type” covers more than blockwise proximal minimization: it also includes stringwise compositions, block averages, and multi-set reflections, provided the fixed-point sets align with the original feasibility problem (Censor et al., 2015).

5. Hadamard-manifold extensions

On symmetric Hadamard manifolds, parallel Douglas–Rachford methods have been developed for ROF-like variational models of the form

$0 \in (A_1+\cdots+A_m)(x)$19

where $0 \in (A_1+\cdots+A_m)(x)$20 and $0 \in (A_1+\cdots+A_m)(x)$21 is a discrete anisotropic total-variation-like regularizer. For a proper, convex, lower-semicontinuous function $0 \in (A_1+\cdots+A_m)(x)$22 on a Hadamard manifold $0 \in (A_1+\cdots+A_m)(x)$23, the proximal map is

$0 \in (A_1+\cdots+A_m)(x)$24

and the reflection is the geodesic symmetry about the proximal point. The Euclidean fact that reflections of proper convex lsc functions are nonexpansive does not carry over in general; the manifold analysis therefore proves nonexpansiveness only for the specific functions needed in the model, namely certain distance-like functions and, for indicator functions, convex sets on constant-curvature symmetric Hadamard manifolds. Under these nonexpansiveness assumptions, the parallel DR update is a Krasnoselski–Mann iteration on $0 \in (A_1+\cdots+A_m)(x)$25: $0 \in (A_1+\cdots+A_m)(x)$26 with $0 \in (A_1+\cdots+A_m)(x)$27, and the projection onto the diagonal yields a minimizer of $0 \in (A_1+\cdots+A_m)(x)$28 (Bergmann et al., 2015).

A later Hadamard-manifold development considers the composite problem

$0 \in (A_1+\cdots+A_m)(x)$29

with each $0 \in (A_1+\cdots+A_m)(x)$30 proper, lower semicontinuous, and geodesically convex. After lifting to $0 \in (A_1+\cdots+A_m)(x)$31 and introducing the diagonal $0 \in (A_1+\cdots+A_m)(x)$32, it defines $0 \in (A_1+\cdots+A_m)(x)$33 and studies three parallel schemes: a non-inertial Krasnoselski–Mann iteration,

$0 \in (A_1+\cdots+A_m)(x)$34

an inertial method,

$0 \in (A_1+\cdots+A_m)(x)$35

and a $0 \in (A_1+\cdots+A_m)(x)$36-accelerated normal $0 \in (A_1+\cdots+A_m)(x)$37-iteration,

$0 \in (A_1+\cdots+A_m)(x)$38

Under the stated parameter conditions, the inertial and $0 \in (A_1+\cdots+A_m)(x)$39-accelerated schemes converge strongly to a fixed point $0 \in (A_1+\cdots+A_m)(x)$40, with residual estimates of order $0 \in (A_1+\cdots+A_m)(x)$41 or $0 \in (A_1+\cdots+A_m)(x)$42 (Sahu et al., 28 Sep 2025).

The Hadamard-manifold literature represented here is not uniform on the scope of reflection nonexpansiveness. One account explicitly states that on general Hadamard manifolds proper convex lower semicontinuous functions can have expansive reflections, whereas the later development formulates its algorithms under the statement that $0 \in (A_1+\cdots+A_m)(x)$43 is nonexpansive on any Hadamard manifold (Bergmann et al., 2015, Sahu et al., 28 Sep 2025). Accordingly, convergence statements on manifolds are tightly linked to the exact nonexpansiveness assumptions adopted in the individual work.

6. Computational structure, applications, and limitations

Across the literature, the computational signature of parallel Douglas–Rachford type algorithms is stable: local resolvent or proximal evaluations are fully decoupled, while global coupling enters only through a diagonal projection, an average, or a primal–dual linear step. In the multioperator Hilbert-space algorithm this means $0 \in (A_1+\cdots+A_m)(x)$44 parallel resolvent evaluations, one all-reduce average, and $0 \in (A_1+\cdots+A_m)(x)$45 local updates (Alcantara et al., 6 Jan 2025). In the primal–dual splitting algorithms it means separate processing of bounded linear operators and set-valued operators at each iteration, with all $0 \in (A_1+\cdots+A_m)(x)$46-indexed computations performed in parallel (Bot et al., 2012). In string- and block-based feasibility schemes, the tradeoff is explicit: longer strings mimic cyclic DR and provide less parallelism, whereas shorter strings or smaller blocks yield more concurrent tasks but slower global progress; that work also notes that no explicit non-asymptotic complexity or rate bounds are given (Censor et al., 2015).

The empirical record in the cited works is correspondingly diverse. For three generalized Heron problems, the two primal–dual algorithms DR1 and DR2 reached 6-digit accuracy in 20 or 50 iterations, with the same final objective values $0 \in (A_1+\cdots+A_m)(x)$47, while subgradient methods require $0 \in (A_1+\cdots+A_m)(x)$48 iterations for comparable accuracy. For image deblurring with

$0 \in (A_1+\cdots+A_m)(x)$49

DR1 and DR2 reduce the objective more rapidly and achieve higher Improvement in Signal-to-Noise Ratio than the forward–backward–forward method of Combettes–Pesquet up to 200 iterations (Bot et al., 2012).

In manifold-valued imaging, the parallel DR algorithm on symmetric Hadamard manifolds reached the lowest TV value, about $0 \in (A_1+\cdots+A_m)(x)$50, in about 300 iterations and about $0 \in (A_1+\cdots+A_m)(x)$51 on a $0 \in (A_1+\cdots+A_m)(x)$52 retinal patch; the cyclic proximal point algorithm required 1500 iterations and about $0 \in (A_1+\cdots+A_m)(x)$53, while half-quadratic minimization reached only a smoothed energy. On SPD$0 \in (A_1+\cdots+A_m)(x)$54 DT-MRI denoising, PDRA converged in about 50 iterations and about $0 \in (A_1+\cdots+A_m)(x)$55, compared with about 1040 iterations and about $0 \in (A_1+\cdots+A_m)(x)$56 for CPPA; HQMA was slightly faster per iteration but did not minimize the true TV energy. In inpainting-plus-denoising on SPD$0 \in (A_1+\cdots+A_m)(x)$57, PDRA completed in 117 iterations and about $0 \in (A_1+\cdots+A_m)(x)$58, compared with 2210 iterations and about $0 \in (A_1+\cdots+A_m)(x)$59 for CPPA (Bergmann et al., 2015). In the later Hadamard-manifold study, the $0 \in (A_1+\cdots+A_m)(x)$60-accelerated scheme outperformed both classic and inertial DR in every generalized Heron test reported; for the Rosenbrock example the iteration counts were 67 for Alg(DR), 32 for Alg(InDR), and 16 for Alg(p-AccDR), and for the Heron problem with $0 \in (A_1+\cdots+A_m)(x)$61 they were 113, 83, and 34, respectively (Sahu et al., 28 Sep 2025).

Several limitations recur. In the manifold-valued ROF setting, rigorous convergence on nonconstant-curvature spaces such as SPD$0 \in (A_1+\cdots+A_m)(x)$62 remains open because nonexpansiveness of the reflection at the diagonal constraint is not guaranteed, even though the numerical tests showed stable convergence (Bergmann et al., 2015). In primal–dual Hilbert-space splitting, the main algorithmic compromise is between larger admissible primal–dual steps and a lower number of evaluations of $0 \in (A_1+\cdots+A_m)(x)$63 and $0 \in (A_1+\cdots+A_m)(x)$64: Algorithm 1 allows $0 \in (A_1+\cdots+A_m)(x)$65 but uses two forward evaluations, whereas Algorithm 2 reduces matrix-vector products and requires $0 \in (A_1+\cdots+A_m)(x)$66 plus an auxiliary variable (Bot et al., 2012). A plausible implication is that the term “parallel Douglas–Rachford type” denotes a family unified less by a single iteration formula than by a common decomposition principle: blockwise resolvents, a consensus mechanism, and convergence analysis through monotone-operator or fixed-point theory.