---
title: Parallel Douglas–Rachford Algorithms
url: https://www.emergentmind.com/topics/parallel-douglas-rachford-type-algorithms
type: topic
---

# Parallel Douglas–Rachford Algorithms

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 [1212.0326, 2501.02752, 1908.05406, 1512.00409, 1512.02814, 2509.23939].

## 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 [2501.02752].

For convex minimization of finitely many functions,
\[
\min_{x\in X}\sum_{i=1}^m g_i(x),
\]
the same principle appears in the product space \(\mathbf X=X^m\), with
\[
g(\mathbf x)=\sum_{i=1}^m g_i(x_i),
\qquad
A=\{(x,\dots,x):x\in X\},
\]
so that the problem becomes
\[
\min_{\mathbf x\in\mathbf X} I_A(\mathbf x)+g(\mathbf x).
\]
The diagonal condition is therefore encoded as a single linear subspace constraint, and the prox of \(g\) splits into the \(m\) independent prox operators \(\prox_{g_i}\) [1908.05406].

On manifolds the same architecture persists. For \(\min_{x\in M}\sum_{k=1}^K \phi_k(x)\) on a symmetric Hadamard manifold, one lifts to \(M^K\) and adds the diagonal set \(\Delta=\{X:x_1=\cdots=x_K\}\). In Hadamard-manifold formulations, the diagonal projection is expressed through a Riemannian centroid,
\[
\bar x=\arg\min_{y\in M}\sum_{i=1}^m d^2(x_i,y),
\]
and \(P_D(x)=(\bar x,\dots,\bar x)\) [1512.02814, 2509.23939].

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 [2501.02752, 1908.05406].

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

For maximal \(\sigma_i\)-monotone operators \(A_i:H\rightrightarrows H\), the block-parallel Douglas–Rachford algorithm on \(H^m\) takes the standard two-operator form
\[
X^{k+1}
=
X^k+\gamma\bigl(J_{\tau G}(2J_{\tau F}(X^k)-X^k)-J_{\tau F}(X^k)\bigr),
\qquad
\gamma\in(0,2).
\]
Because \(F\) and \(G\) are separable, the iteration can be written as three explicit steps:
\[
U_i^k=J_{\tau A_i}(X_i^k)\quad(i=1,\dots,m)\ \text{in parallel},
\]
\[
V_i^k=\bar y^k,
\qquad
\bar y^k=\frac1m\sum_{j=1}^m(2U_j^k-X_j^k),
\]
\[
X_i^{k+1}=X_i^k+\gamma(\bar y^k-U_i^k)\quad(i=1,\dots,m)\ \text{in parallel}.
\]
The “shadow” sequence \(Z^k=U^k\) and the projected sequence \(Y^k=V^k\) both serve as candidate solution sequences [2501.02752].

The convergence theory is divided into regimes by the aggregate monotonicity parameter \(\Sigma=\sum_i \sigma_i\). In the convex/monotone case, where all \(\sigma_i=0\), any \(\tau>0\) and any \(\gamma\in(0,2)\) yield the classical Douglas–Rachford behavior: \(X^k\rightharpoonup X^*\in\mathrm{Fix}(R)\), and \(Z^k\to(x^*,\dots,x^*)\) with \(x^*\in\operatorname{zer}(\sum A_i)\). If \(\Sigma>0\), then one chooses \(\tau\) small enough so that \(1+\tau\min_i\{\sigma_i\}>0\); in that case the method converges strongly, \(x^*\) is unique, and both \(X^k\) and \(Z^k\) converge to it. In the finite-dimensional nonconvex setting considered in the same work, where \(m-1\) operators are Lipschitz gradients of \(C^1\) functions and the remaining one is a proper closed subdifferential, blockwise choices \(\tau<\min_i(1/L_i)\) and \(\gamma\in(0,2)\) guarantee that every cluster point of the shadow sequence is a stationary point, with residual \(\|(I-R)X^k\|=o(1/\sqrt{k})\) [2501.02752].

The implementation pattern is explicit: each sweep consists of \(m\) parallel resolvent calls, one global sum/average, and \(m\) 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 [2501.02752].

## 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 \(\mathcal H\) and \(\mathcal H_i\), with maximally monotone operators \(A:\mathcal H\to2^\mathcal H\), \(B_i,D_i:\mathcal H_i\to2^{\mathcal H_i}\), bounded linear operators \(L_i:\mathcal H\to\mathcal H_i\), and data \(z\in\mathcal H\), \(r_i\in\mathcal H_i\), the primal inclusion is formulated as
\[
0 \in A x + \sum_{i=1}^m L_i^*(B_i\, D_i)(L_i x-r_i),
\]
and equivalently, writing \(M_i\) for the parallel sum of \(B_i\) and \(D_i\),
\[
z \in A x + \sum_{i=1}^m L_i^* M_i(L_i x-r_i).
\]
The associated dual inclusion is
\[
-z+\sum_{i=1}^m L_i^*v_i\in A x,
\qquad
v_i\in(B_i\, D_i)(L_i x-r_i),\quad i=1,\dots,m.
\]
Under the qualification
\[
z\in\operatorname{Ran}\Bigl(A+\sum_i L_i^*(B_i\, D_i)(L_i\bullet-r_i)\Bigr),
\]
the set of primal–dual solutions is nonempty [1212.0326].

The paper develops two inexact Douglas–Rachford type primal–dual algorithms. Algorithm 1 is a two-step Douglas–Rachford splitting obtained by lifting \((A,S,Q)\) to a product space \(\mathcal K\) and applying an inexact DR iteration in a metric induced by a strongly positive self-adjoint operator \(V\). With step sizes \(\tau>0\) and \(\sigma_i>0\) satisfying
\[
\tau\sum_{i=1}^m \sigma_i \|L_i\|^2 < 4,
\]
relaxation parameters \(\lambda_n\in(0,2)\), and summable error sequences, the method uses the resolvents \(J_{\tau A}\), \(J_{\sigma_i B_i^{-1}}\), and \(J_{\sigma_i D_i^{-1}}\); all computations for \(i=1,\dots,m\) can be carried out in parallel. Algorithm 2 is a single-step Douglas–Rachford splitting in which each evaluation of \(L_i\) and \(L_i^*\) is reduced to a single forward–backward pass. Its price is a stricter condition,
\[
\tau\sum_{i=1}^m \sigma_i \|L_i\|^2 < \frac14,
\]
together with an auxiliary variable \(y_{i,n}\) [1212.0326].

The convergence results mirror the operator-theoretic structure. For Algorithm 1, if \(A\) and \(B_i^{-1}\) are merely maximally monotone, then \((x_n,v_{1,n},\dots,v_{m,n})\) converges weakly to a primal–dual solution. In finite dimensions, if \(\lambda_n\) is bounded away from \(0\) and \(A,B_i^{-1}\) are uniformly monotone, then \(x_n\to x^*\) and \(v_{i,n}\to v_i^*\) strongly. For Algorithm 2, weak convergence holds under the same range condition, while strong convergence in finite dimensions requires uniform monotonicity of \(A\), \(B_i^{-1}\), and \(D_i\) [1212.0326].

The practical comparison between the two schemes is explicit. Algorithm 1 uses two forward evaluations of \(L_i\) and \(L_i^*\) per iteration but allows the larger stepsize bound \(\tau\sum_i\sigma_i\|L_i\|^2<4\). Algorithm 2 halves the matrix-vector products at the cost of the stricter bound \(\tau\sum_i\sigma_i\|L_i\|^2<\frac14\) and the additional auxiliary variable \(y_{i,n}\). Both treat the resolvents in parallel and admit inexact computations via summable errors [1212.0326].

## 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
\[
\min_{x\in X}\sum_{i=1}^m g_i(x),
\]
the product-space reformulation with diagonal subspace \(A\subset X^m\) yields the parallel componentwise update
\[
x_{n+1,i}
=
x_{n,i}-\bar x_n+\prox_{g_i}(2\bar x_n-x_{n,i}),
\qquad
\bar x_n=\frac1m\sum_{i=1}^m x_{n,i}.
\]
This is precisely the Douglas–Rachford operator \(T=\Id-P_A+\prox_g\circ R_A\) written in coordinates, so each iteration consists of \(m\) 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^\perp+\dom g^*\), and nonemptiness of the normal solution set, the shadow sequence \(\bar x_n=P_A T^n(\mathbf x_0)\) converges to a limit \(\bar x\) solving
\[
\min_x \sum_{i=1}^m g_i(x-v_i),
\]
where \(v=(v_1,\dots,v_m)=P_{A^\perp-\dom g}(0)\). This is the mechanism that allows the method to remain meaningful even in possibly inconsistent cases [1908.05406].

For convex feasibility with closed convex sets \(C_i\subset\mathcal H\), new algorithmic structures embed the basic two-set Douglas–Rachford operator \(T_{B,A}=\frac12(\Id+R_A R_B)\) into larger parallel architectures. In String-Averaging DR, one partitions the index set into strings \(I_t\), composes two-set DR operators along each string, and averages the resulting string operators with weights \(w_t\). If \(\operatorname{int}(\bigcap_i C_i)\neq\emptyset\), the sequence converges strongly to a point in \(\bigcap_i C_i\). 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 \(\bigcap_i C_i\neq\emptyset\), the swept iterates converge weakly, and strong convergence follows under the interior-point condition. The same work also defines an \(r\)-set DR operator
\[
T_{C_1,\dots,C_r}=\tfrac12(\Id+V_{C_1,\dots,C_r}),
\]
where \(V_{C_1,\dots,C_r}\) is a composite reflection, and studies mixtures
\[
x^{k+1}=\sum_{r=2}^m w_r\,T_{C_1,\dots,C_r}(x^k).
\]
These schemes include the cyclic Douglas–Rachford algorithm, the simultaneous DR operator, and averaged DR as special cases [1512.00409].

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 [1512.00409].

## 5. Hadamard-manifold extensions

On symmetric Hadamard manifolds, parallel Douglas–Rachford methods have been developed for ROF-like variational models of the form
\[
E(u)=D(u;f)+\alpha\,TV(u),
\]
where \(D(u;f)=\frac12\sum_i d(u_i,f_i)^2\) and \(TV(u)\) is a discrete anisotropic total-variation-like regularizer. For a proper, convex, lower-semicontinuous function \(\phi\) on a Hadamard manifold \(M\), the proximal map is
\[
\Prox_{\lambda\phi}(x)=\arg\min_{y\in M}\Bigl\{\phi(y)+\frac1{2\lambda}d(x,y)^2\Bigr\},
\]
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 \(M^K\):
\[
X^{(r+1)}=\gamma_r S^{(r)}\oplus (1-\gamma_r)X^{(r)},
\qquad
S^{(r)}=R_{\eta\Phi}\circ R_{\eta I_\Delta}(X^{(r)}),
\]
with \(\sum_r \gamma_r(1-\gamma_r)=\infty\), and the projection onto the diagonal yields a minimizer of \(\sum_k \phi_k\) [1512.02814].

A later Hadamard-manifold development considers the composite problem
\[
\min_{x\in M} F(x),\qquad F(x)=\sum_{i=1}^m f_i(x),
\]
with each \(f_i\) proper, lower semicontinuous, and geodesically convex. After lifting to \(M^m\) and introducing the diagonal \(D\), it defines \(T=R_{\gamma\widehat F}\circ R_{\gamma\iota_D}\) and studies three parallel schemes: a non-inertial Krasnoselski–Mann iteration,
\[
x^{k+1}=\gamma(x^k,Tx^k;\alpha_k),
\]
an inertial method,
\[
y^k=\exp_{x^k}\!\bigl(-\theta_k\exp_{x^k}^{-1}(x^{k-1})\bigr),\qquad
x^{k+1}=\gamma(y^k,Ty^k;\alpha_k),
\]
and a \(p\)-accelerated normal \(S\)-iteration,
\[
y^k=\gamma(x^k,Tx^k;\alpha_k),\qquad x^{k+1}=T^p(y^k).
\]
Under the stated parameter conditions, the inertial and \(p\)-accelerated schemes converge strongly to a fixed point \(v\in\Fix T\), with residual estimates of order \(O(1/\sqrt{k})\) or \(O(1/\sqrt{\sum_{i=1}^k\alpha_i(1-\alpha_i)})\) [2509.23939].

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 \(R_{\gamma f}\) is nonexpansive on any Hadamard manifold [1512.02814, 2509.23939]. 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 \(m\) parallel resolvent evaluations, one all-reduce average, and \(m\) local updates [2501.02752]. In the primal–dual splitting algorithms it means separate processing of bounded linear operators and set-valued operators at each iteration, with all \(i\)-indexed computations performed in parallel [1212.0326]. 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 [1512.00409].

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 \(V^*\), while subgradient methods require \(10^6+\) iterations for comparable accuracy. For image deblurring with
\[
\min_{x\in[0,1]^n}\|Ax-b\|_1+\alpha_2\|Wx\|_1+\alpha_1 TV(x),
\]
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 [1212.0326].

In manifold-valued imaging, the parallel DR algorithm on symmetric Hadamard manifolds reached the lowest TV value, about \(184.36\), in about 300 iterations and about \(34\,\mathrm{s}\) on a \(50\times50\) retinal patch; the cyclic proximal point algorithm required 1500 iterations and about \(57\,\mathrm{s}\), while half-quadratic minimization reached only a smoothed energy. On SPD\((3)\) DT-MRI denoising, PDRA converged in about 50 iterations and about \(140\,\mathrm{s}\), compared with about 1040 iterations and about \(640\,\mathrm{s}\) for CPPA; HQMA was slightly faster per iteration but did not minimize the true TV energy. In inpainting-plus-denoising on SPD\((3)\), PDRA completed in 117 iterations and about \(36\,\mathrm{s}\), compared with 2210 iterations and about \(161\,\mathrm{s}\) for CPPA [1512.02814]. In the later Hadamard-manifold study, the \(p\)-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 \(m=2,N=4\) they were 113, 83, and 34, respectively [2509.23939].

Several limitations recur. In the manifold-valued ROF setting, rigorous convergence on nonconstant-curvature spaces such as SPD\((3)\) remains open because nonexpansiveness of the reflection at the diagonal constraint is not guaranteed, even though the numerical tests showed stable convergence [1512.02814]. In primal–dual Hilbert-space splitting, the main algorithmic compromise is between larger admissible primal–dual steps and a lower number of evaluations of \(L_i\) and \(L_i^*\): Algorithm 1 allows \(\tau\sum_i \sigma_i\|L_i\|^2<4\) but uses two forward evaluations, whereas Algorithm 2 reduces matrix-vector products and requires \(\tau\sum_i \sigma_i\|L_i\|^2<\frac14\) plus an auxiliary variable [1212.0326]. 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.

Source: https://www.emergentmind.com/topics/parallel-douglas-rachford-type-algorithms