---
title: Distributed In-Memory PDHG Techniques
url: https://www.emergentmind.com/topics/distributed-in-memory-primal-dual-hybrid-gradient-pdhg
type: topic
---

# Distributed In-Memory PDHG Techniques

Distributed in-memory Primal–Dual Hybrid Gradient (PDHG) denotes a class of large-scale first-order optimization methods and implementations built around repeated applications of a linear operator and its transpose, together with proximal or projection steps. In the recent literature, the term spans several closely related developments: exact primal-only reformulations of PDHG for linearly constrained convex programs, matrix-free PDHG solvers for large sparse linear programming, distributed multi-GPU realizations in which the matrix and working vectors remain resident in aggregate GPU memory, and RRAM-based in-memory implementations that encode a fixed operator once and reuse it throughout the solve [1706.02602] [2106.04756] [2601.07628] [2509.21137].

## 1. Problem classes and saddle-point formulations

A foundational formulation is the linearly constrained convex program
$$
\min_{x\in \mathbb{R}^n} \quad g(x) \quad \text{s.t.} \quad Ax=b,
$$
where \(A\in \mathbb{R}^{m\times n}\), \(b\in \mathbb{R}^m\), and \(g:\mathbb{R}^n\to(-\infty,+\infty]\) is proper, lower semicontinuous, convex. Introducing a dual multiplier \(y\in\mathbb{R}^m\), the problem is written as
$$
\min_{x\in \mathbb{R}^n}\max_{y\in\mathbb{R}^m}\quad g(x)+\langle Ax,y\rangle-\langle b,y\rangle.
$$
For this class, PDHG can be written as an entirely primal algorithm, and the reformulation remains meaningful even in degenerate cases when the linear system is inconsistent or when strong duality does not hold [1706.02602].

A central auxiliary object is
$$
f(x)=\frac12\|Ax-b\|^2,\qquad f_*:=\min_x f(x),
$$
which induces the generalized constrained problem
$$
\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad x\in \arg\min f.
$$
Equivalently,
$$
\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad A^\top A x=A^\top b.
$$
This covers consistent systems, inconsistent systems, and cases where strong duality may fail. The significance for distributed implementations is direct: real systems often involve noisy, ill-conditioned, or only approximately consistent constraints, so convergence statements that do not rely on exact saddle-point existence enlarge the operational scope of PDHG [1706.02602].

In linear programming, the same PDHG logic appears in bounded-form and box-constrained saddle formulations. One distributed PDLP line considers
$$
\min_{x \in \mathcal{X}} c^\top x \quad \text{subject to } Ax \in \mathcal{S},
$$
with
$$
\mathcal{X} := \{x \in \mathbb{R}^n : \ell_v \le x \le u_v\}, \qquad \mathcal{S} := \{s \in \mathbb{R}^m : \ell_c \le s \le u_c\},
$$
and the saddle function
$$
\min_{x \in \mathcal{X}} \max_{y \in \mathcal{Y}} \; \mathcal{L}(x,y) \coloneqq c^\top x + y^\top A x - p(-y; \ell_c, u_c).
$$
In this LP specialization, the dominant kernels remain sparse matrix-vector products \(Ax\) and \(A^\top y\), together with projections onto box sets, which is precisely the computational anatomy that motivates distributed-memory and in-memory deployments [2601.07628].

## 2. Exact primal-only reduction and decentralized consensus

For the saddle problem
$$
\min_x \max_y \; g(x)+\langle Ax,y\rangle-\langle b,y\rangle,
$$
the specialized PDHG iteration is
$$
\begin{cases}
y^{k+1}=y^k+\sigma(A\bar x^k-b),\\[2mm]
x^{k+1}=\prox_{\tau g}\bigl(x^k-\tau A^\top y^{k+1}\bigr),
\end{cases}
\qquad \bar x^k=2x^k-x^{k-1}.
$$
Under \(y^0=0\) and \(\bar x^0=x^0\), eliminating the dual recursion yields the exact primal-only form
$$
\begin{cases}
z^k=\dfrac{x^k+ks^k}{k+1},\\[2mm]
x^{k+1}=\prox_{\tau g}\bigl(x^k-(k+1)\lambda\nabla f(z^k)\bigr),\\[2mm]
s^{k+1}=\dfrac{x^{k+1}+ks^k}{k+1},
\end{cases}
\qquad f(x)=\frac12\|Ax-b\|^2,\quad \lambda=\tau\sigma,
$$
and, using linearity of \(\nabla f\),
$$
\begin{cases}
x^{k+1}=\prox_{\tau g}\bigl(x^k-\lambda \nabla f(x^k+ks^k)\bigr),\\[2mm]
s^{k+1}=\dfrac{x^{k+1}+ks^k}{k+1}.
\end{cases}
$$
The sequence \((x^k)\) generated by the primal-only method is identical to that produced by PDHG, under the simplifying initialization \(y^0=0\). The method is interpreted as a modified Tseng accelerated proximal-gradient scheme with a penalty parameter increasing to infinity; with
$$
F_k(x)=g(x)+k\lambda\bigl(f(x)-f_*\bigr),
$$
it can be viewed as performing one proximal-gradient-like step for each increasing-penalty problem \(\min_x F_k(x)\) [1706.02602].

The convergence theory is broader than standard saddle-point analyses. Under \(S\neq\varnothing\), \(\tau,\sigma>0\), and \(\tau\sigma\|A\|^2<1\), the averaged sequence satisfies
$$
F_k(s^k)-g_*=O(1/k),
$$
more precisely
$$
F_k(s^k)-g_*\le \frac{D_x^2}{2k}.
$$
If strong duality holds, then both \(x^k\) and \(s^k\) converge to a solution of the generalized problem, with
$$
f(x^k)-f_*=O(1/k), \qquad f(s^k)-f_*=O(1/k^2), \qquad |g(s^k)-g_*|=O(1/k).
$$
Without strong duality, but with bounded solution set and \(g\) bounded below, all cluster points of \((s^k)\) belong to \(S\), and
$$
f(s^k)-f_*=o(1/k).
$$
This changes the interpretation of PDHG in degenerate regimes: the averaged primal iterates continue to solve a meaningful generalized problem even when the original feasibility system is inconsistent [1706.02602].

The decentralized specialization is particularly important for distributed in-memory settings. For consensus over a connected undirected graph \(G=(V,E)\),
$$
\min_{x_1,\dots,x_n} g_1(x_1)+\cdots+g_n(x_n) \quad \text{s.t.}\quad x_1=\cdots=x_n,
$$
stacking \(x=(x_1,\dots,x_n)\) gives
$$
\min_x g(x)\quad \text{s.t.}\quad Lx=0,
$$
where \(L\) is the graph Laplacian. Replacing the equivalent constraint \(Lx=0\) with \(\sqrt L\,x=0\) and then applying the primal-only reduction yields
$$
\begin{cases}
x^{k+1}=\prox_{\tau g}\bigl(x^k-\lambda L(x^k+ks^k)\bigr),\\[1mm]
s^{k+1}=\dfrac{x^{k+1}+ks^k}{k+1}.
\end{cases}
$$
Per node,
$$
x_i^{k+1} = \prox_{\tau g_i}\Bigl( x_i^k-\lambda \bigl[L(x^k+ks^k)\bigr]_i \Bigr),
$$
with a distributed implementation:
1. form \(w_i^k=x_i^k+ks_i^k\);
2. exchange \(w_i^k\) with neighbors;
3. compute \(\ell_i^k=\deg(v_i)w_i^k-\sum_{j\in\mathcal N_i} w_j^k\);
4. update \(x_i^{k+1}=\prox_{\tau g_i}(x_i^k-\lambda \ell_i^k)\);
5. update \(s_i^{k+1}=\frac{x_i^{k+1}+ks_i^k}{k+1}\).

Only one neighborhood communication is needed. By contrast, direct PDHG on \(Lx=0\) requires \(L\bar x^k\) and \(Ly^{k+1}\), hence two communication rounds per iteration, and storage of \(x_i^k\), \(x_i^{k-1}\), and \(y_i^k\). The primal-only scheme stores only \(x_i^k\) and \(s_i^k\), and the step-size condition improves from \(\tau\sigma \|L\|^2<1\) to \(\tau\sigma \|L\|<1\) [1706.02602].

## 3. PDLP and matrix-free large-scale linear programming

A second major strand is PDLP, a practical first-order method for linear programming derived by applying PDHG to a saddle-point formulation of LP. The primal-dual pair considered is
$$
\begin{aligned}
\min_{x \in \mathbb{R}^n}\quad & c^\top x \\
\text{s.t.}\quad & Gx \ge h, \\
& Ax = b, \\
& l \le x \le u,
\end{aligned}
$$
with the associated saddle problem
$$
\min_{x \in X}\max_{y \in Y} \; c^\top x - y^\top Kx + q^\top y,
$$
where
$$
X := \{x \in \mathbb{R}^n : l \le x \le u\}, \qquad Y := \{y \in \mathbb{R}^{m_1+m_2} : y_{1:m_1} \ge 0\}.
$$
The baseline PDHG iteration is
$$
\begin{aligned}
x^{k+1} &= \Pi_X\!\left(x^k - \tau \left(c - K^\top y^k\right)\right), \\
y^{k+1} &= \Pi_Y\!\left(y^k + \sigma \left(q - K(2x^{k+1}-x^k)\right)\right),
\end{aligned}
$$
with the reparameterization
$$
\tau = \frac{\eta}{\omega}, \qquad \sigma = \omega \eta,
$$
so the standard convergence condition becomes \(\eta \le \frac{1}{\|K\|_2}\) [2106.04756].

PDLP is not just one fixed PDHG iteration; it is an outer/inner restarted scheme around PDHG with adaptive step selection and primal-weight updates, preceded by presolve and diagonal scaling. The enhancements include diagonal preconditioning, presolving, adaptive step sizes, adaptive restarting, and primal/dual weighting. The method is explicitly matrix-free in the sense that it never factorizes the KKT system. Each inner iteration performs one \(K^\top y\), one \(K(\cdot)\) on an extrapolated primal vector, simple projections, and scalar or norm computations for step-size acceptance. For large sparse LPs, the core operation is therefore repeated sparse matrix-vector multiplication, with low auxiliary state and no factorization storage [2106.04756].

The diagonal scaling pipeline is operationally important. PDLP uses 10 iterations of Ruiz scaling, then Pock–Chambolle scaling. This balances magnitudes of rows and columns, reduces effective conditioning for first-order updates, stabilizes a single global scalar step size, and makes primal and dual progress more symmetric. Presolve transformations include detecting inconsistent bounds, removing empty rows or columns, removing fixed variables, detecting duplicate rows, and tightening bounds. In a distributed-memory context, these choices are attractive because they preserve sparsity pattern while reducing memory footprint and matvec cost [2106.04756].

The convergence theory is stratified. For the LP saddle problem, baseline PDHG converges to an optimal solution if
$$
\tau \sigma \|K\|_2^2 \le 1.
$$
Restart theory is partially grounded, and the paper states that sufficient normalized-gap decay guarantees the linear convergence of restarted PDHG on LP problems under the assumptions of the cited restart work. At the same time, PDLP as a whole has no convergence proof: adaptive step size has no proof of convergence in the paper, primal weight updates are not proved, and the final production combination is partly heuristic [2106.04756].

Empirically, PDLP is positioned as a high-accuracy first-order LP solver. On a set of 383 LP instances derived from MIPLIB 2017, with a target of \(10^{-8}\) relative accuracy and 1 hour time limit, PDLP achieves a 6.3x reduction in the geometric mean of solve times and a 4.6x reduction in the number of instances unsolved, from 227 to 49, relative to SCS. On the MIP Relaxations benchmark at \(\epsilon=10^{-8}\), PDLP solved 334 instances with shifted geometric mean time 106.4s, compared with 80 for baseline PDHG and 124 for matrix-free SCS. The PageRank study reports 7.4s at \(10^4\) nodes, 35s at \(10^5\) nodes, 11 min at \(10^6\) nodes, and 5.4 hr at \(10^7\) nodes; at \(10^6\) nodes, Gurobi barrier is out of memory. These results do not by themselves define a distributed system, but they establish the matrix-free, SpMV-dominated operating regime that later distributed-memory PDHG implementations exploit [2106.04756].

## 4. Step-size control, convergence monitoring, and limits of ratio tuning

Step-size choice is a persistent issue in PDHG because the admissible region, such as \(\sigma\tau\|A\|^2<1\), is usually broad while asymptotic speed can vary substantially inside it. One recent line of work studies adaptive selection through convergence-rate monitoring rather than only through residual balancing. The problem class is
$$
\min_{x\in\mathbb R^n} f(x) + f_2(x) + g(Ax),
$$
rewritten as the saddle-point problem
$$
\min_{x\in\mathbb R^n}\max_{y\in\mathbb R^m} L(x,y) = f(x)+f_2(x)+\langle Ax,y\rangle-g^*(y).
$$
For the Vũ–Condat form of PDHG,
$$
y_{k+1} = \prox_{\sigma_k g^*}\big(y_k + \sigma_k A x_k\big), \qquad
x_{k+1} = \prox_{\tau_k f}\Big(x_k - \tau_k \nabla f_2(x_k) - \tau_k A^\top(2y_{k+1}-y_k)\Big),
$$
with
$$
\sigma_k \tau_k \|A\|^2 + \tau_k L_f/2 < 1,
$$
the paper shows how to estimate the spectral radius of the linearized iteration from residual norms alone and to tune \((\tau,\sigma)\) in order to reduce that spectral radius [2403.19202].

The central estimator is based on
$$
r_k := \frac{\|z_{k+1}-z_k\|_{V_s}}{\|z_k-z_{k-1}\|_{V_s}},
$$
together with a procedure that remains effective even when the dominant eigenvalues are a complex conjugate pair. The method uses only consecutive iterate differences, the metric norm \(\|\cdot\|_{V_s}\), and local checks on the scalar sequence \(r_k\). This is technically important for distributed in-memory PDHG because it replaces explicit operator-level spectral analysis with norm-based global diagnostics. For strongly convex quadratics, the proposed step-size rule yields an algorithm as fast as inertial gradient descent [2403.19202].

Residual balancing remains relevant, especially as a warm-up phase. Goldstein’s rule compares primal and dual residual magnitudes and shifts the ratio \((\tau,\sigma)\) while preserving roughly the product. The convergence-monitoring paper’s conclusion is more nuanced: residual balance reacts quickly early on and detects gross primal/dual imbalance well, but it often stalls before reaching the truly best asymptotic ratio, especially when iterates oscillate due to complex dominant eigenvalues. This suggests a layered policy for large-scale implementations: use residual balancing to enter a reasonable regime, then use rate monitoring for asymptotic tuning [2403.19202].

The same paper also imposes an important limitation. For linear programs, asymptotic PDHG rate depends only on the product \(\sigma\tau\), not their ratio, so ratio adaptation may not improve asymptotic speed there. The reported sparse SVM or LP-like experiments show that changing the primal/dual ratio does not improve asymptotic rate, and that averaging and restarting PDHG provides more gain than step-size ratio tuning. In distributed LP solvers, this narrows the expected benefit of sophisticated ratio adaptation: it can still affect the transient active-set identification phase, but it is not a universal asymptotic accelerator [2403.19202].

## 5. Distributed multi-GPU PDLP

A direct distributed-memory realization appears in work on scaling PDLP beyond a single GPU. The target is massive-scale LPs solved by a distributed-memory, multi-GPU implementation of PDHG/PDLP built on top of cuPDLPx. The motivating setting is LPs with millions to hundreds of millions of variables or constraints and up to hundreds of millions of nonzeros, where a single GPU is limited by both memory capacity and throughput. The distributed implementation aggregates memory across GPUs so the matrix and vectors can remain resident in GPU memory, and it aggregates compute and memory bandwidth for the SpMV-dominated iteration [2601.07628].

The central system idea is a 2D process or device grid with row axis \(R\) and column axis \(C\). The sparse matrix \(A\in\mathbb{R}^{m\times n}\) is partitioned into blocks \(A_{[i,j]}\). Vectors in \(\mathbb{R}^n\), including \(x\), \(c\), and variable bounds, are partitioned by the column axis and replicated vertically across all devices in a given column. Vectors in \(\mathbb{R}^m\), including \(y\) and row bounds, are partitioned by the row axis and replicated horizontally across all devices in a given row. On device \((i,j)\), local memory contains \(A_{[i,j]}\), the local primal block \(x_{[j]}\), \(c_{[j]}\), bounds for that column block, and the local dual block \(y_{[i]}\) with row bounds. This lets each device compute \(A_{[i,j]}x_{[j]}\), \(A_{[i,j]}^\top y_{[i]}\), local primal projection, and local dual projection without gathers [2601.07628].

One PDHG iteration has two synchronized vector reductions. First, each device computes
$$
g_{[i,j]} = A_{[i,j]}^\top y_{[i]}^k,
$$
then an AllReduce over \(R\) yields
$$
[A^\top y^k]_{[j]} = \operatorname{AllReduce}_R\bigl(A_{[i,j]}^\top y_{[i]}^k\bigr),
$$
followed by the local primal update
$$
x_{[j]}^{k+1} = \operatorname{proj}_{\mathcal{X}_{[j]}}
\left( x_{[j]}^k - \tau\bigl(c_{[j]} - [A^\top y^k]_{[j]}\bigr) \right).
$$
Second, each device forms \(\bar x = 2x^{k+1}-x^k\), computes
$$
v_{[i,j]} = A_{[i,j]}\bar{x}_{[j]},
$$
and an AllReduce over \(C\) yields
$$
z_{[i]} = \operatorname{AllReduce}_C\bigl(A_{[i,j]}\bar{x}_{[j]}\bigr),
$$
followed by
$$
y_{[i]}^{k+1} = y_{[i]}^k - \sigma z_{[i]} - \sigma \cdot \operatorname{proj}_{-\mathcal{S}_{[i]}}
\left(\sigma^{-1} y_{[i]}^k - z_{[i]}\right).
$$
The reflected Halpern step is fully local because current iterate and anchor are partitioned identically [2601.07628].

| Operation | Collective | Frequency |
|---|---|---|
| Primal step \(A^\top y\) | AllReduce over \(R\) | 1/iter |
| Dual step \(Ax\) | AllReduce over \(C\) | 1/iter |
| Halpern update | none | none |
| KKT evaluation | vector sums and scalar reductions | every \(1/K\) iterations |

This design reduces each communication to a reduction over one axis only, with communicated vector sizes approximately \(O(n/|C|)\) for \(A^\top y\) and \(O(m/|R|)\) for \(Ax\). The implementation uses NCCL AllReduce for GPU-resident vector reductions and MPI AllReduce for small CPU-resident scalar reductions. Load balancing is handled through block-wise random shuffling combined with nonzero-aware data distribution. Natural ordering can produce catastrophic imbalance; on zib03 with 8 GPUs, no permutation and no NNZ balancing required 7368.89s, full random required about 300s, and block random required about 244–245s. The shifted geometric mean over the tested set improved from 220.08s for no permutation and no NNZ balancing to 75.85s for block random plus NNZ balancing [2601.07628].

The reported experiments were run on 1 node with 8 NVIDIA H100 GPUs, each with 80 GB HBM3 and NVLink interconnect, so the measured system is distributed in design but intra-node in evaluation. Scaling is strongly problem-size dependent. For MIPLIB Small, 1 GPU achieved 5.81s while 8 GPUs took 6.80s; for MIPLIB Medium, performance was roughly flat; for MIPLIB Large, time fell from 115.68s on 1 GPU to 66.97s on 8 GPUs, and iteration rate rose from 1303 iter/s to 2758 iter/s. On larger instances, zib03 improved from 812s on 1 GPU to 245s on 8 GPUs; \(sdm\_50k\_500k\_15\_10\) improved from 377s to 61s; and \(mcf\_5000\_100\_250\) improved from 9173s to 1732s. The implementation preserves full FP64 numerical accuracy, but the paper also states that vector synchronization is the main communication bottleneck and that small instances can regress because communication overhead outweighs compute gains [2601.07628].

## 6. RRAM-based distributed in-memory PDHG

A different meaning of distributed in-memory PDHG appears in work on RRAM in-memory computing. Here the concern is not only distributed memory across processors but also co-design for analog crossbar arrays where matrix programming is expensive in time and energy, write endurance is limited, and repeated matrix reprogramming is undesirable. The target problem is a large-scale constrained LP, reduced in implementation to standard form,
$$
\begin{aligned}
\min_{\mathbf{x}\in\mathbb{R}^n}\quad & \mathbf{c}^\top \mathbf{x} \\
\text{s.t.}\quad & \mathbf{K}\mathbf{x} = \mathbf{b},\quad x_i \geq 0.
\end{aligned}
$$
The dominant operations are again matrix–vector multiplications, but the hardware cost model changes completely: the algorithm is redesigned so that the matrix is encoded once, all iterations reuse the same physical mapping, and only vectors are changed iteratively [2509.21137].

The key hardware mapping is the symmetric block matrix
$$
\mathbf{M} = \begin{bmatrix} \mathbf{0}_{m\times m} & \mathbf{K} \\
\mathbf{K}^\top & \mathbf{0}_{n\times n} \end{bmatrix},
$$
which satisfies
$$
\mathbf{M} \begin{bmatrix} \mathbf{y} \\ \mathbf{x} \end{bmatrix}
=
\begin{bmatrix} \mathbf{K}\mathbf{x} \\ \mathbf{K}^\top\mathbf{y} \end{bmatrix}.
$$
This permits one static encoded operator to serve both forward and transpose multiplies, avoiding separate programming of \(\mathbf{K}\) and \(\mathbf{K}^\top\). The accelerator abstraction supports a single device-side MVM kernel with padded inputs for `full`, `A@x`, and `AT@y`. The physical architecture is a 4-by-4 array of 64-by-64 crossbars, acting as a logical 256-by-256 matrix engine; the matrix is partitioned across the 16 crossbars, encoded once using a write-and-verify scheme, and then reused for both Lanczos norm estimation and PDHG iterations [2509.21137].

The implemented algorithm is a custom, hardware-aware enhanced PDHG rather than textbook Chambolle–Pock. Preprocessing performs Ruiz rescaling, diagonal preconditioning, and symmetric-block construction. The operator norm is estimated by Lanczos on \(\mathbf{M}_d\), and initialization uses
$$
\tau = \sigma = \eta / \hat{\sigma}_{\max}(\mathbf{K}), \qquad \eta = 0.95.
$$
Each iteration then uses a deterministic Nesterov-style factor
$$
\theta_k = \frac{1}{\sqrt{1 + 2 \gamma \tau}},
\qquad
\tau \leftarrow \theta_k \tau,
\qquad
\sigma \leftarrow \sigma/\theta_k,
$$
followed by exactly two accelerator MVMs per PDHG iteration: one for \(\mathbf{K}\bar{\mathbf{x}}^{\{k\}}\) and one for \(\mathbf{K}^\top \mathbf{y}^{\{k+1\}}\). The design explicitly rejects line-search-based adaptation because multiple trial MVMs per outer iteration are undesirable on IMC hardware [2509.21137].

Theoretical support is provided for both norm estimation and inexact PDHG under analog perturbations. For
$$
\mathbf{M}:= \begin{bmatrix} \mathbf{0}_{m\times m} & \mathbf{K} \\
\mathbf{K}^\top & \mathbf{0}_{n\times n} \end{bmatrix},
$$
the dominant eigenvalue of \(\mathbf{M}\) equals the dominant singular value of \(\mathbf{K}\),
$$
\lambda_{\max}(\mathbf{M}) \equiv \sigma_{\max}(\mathbf{K}),
$$
which justifies Lanczos on the symmetric block matrix. Under noisy MVMs, the paper gives
$$
\mathbb{E}|\bar\theta_K-L| \le \mathcal{O}\!\left(\frac{1}{K}\right) + \mathcal{O}(\epsilon_{\max}),
$$
and, for the ergodic average of the noisy PDHG iterates,
$$
\mathbb{E}[\mathrm{gap}(\bar{\mathbf{z}}_K)] \le \frac{C_0}{K} + \frac{\delta}{\sqrt{K}}.
$$
The perturbation model assumes unbiased, bounded, finite-variance errors in the effective operator applications. This does not eliminate the asymptotic effect of analog noise, but it preserves the standard convergence structure up to an additive hardware-noise term [2509.21137].

The RRAM evaluation uses MELISO+ as the full-stack simulation platform, with EpiRAM and TaO\(_x\)-HfO\(_x\) device technologies. Relative to GPU PDHG, the reported improvements are largest in the PDHG phase, which dominates total energy and latency. On neos5, PDHG energy savings are 768.12x for EpiRAM and 4999.18x for TaO\(_x\)-HfO\(_x\); PDHG latency speedups are 53.51x and 323.61x, respectively. End-to-end total-process improvements over GPU reach 521x power and 37x latency for EpiRAM on neos5, and 1813x power and 137x latency for TaO\(_x\)-HfO\(_x\) on the same benchmark. Accuracy is described as comparable, but not uniformly at the \(10^{-5}\) level: for gen-ip021, GPU achieved \(2.64\times 10^{-6}\) while both RRAM technologies were approximately \(2.98\times 10^{-2}\), and for gen-ip054 the RRAM gaps were \(3.36\times 10^{-3}\) and \(6.98\times 10^{-3}\). The paper therefore supports both the promise and the limits of analog in-memory PDHG: it demonstrates the first PDHG-based LP solver implemented on RRAMs, but it is simulation-based, benchmark sizes are constrained by the logical array capacity, and some low-level encoding details are omitted in the main text [2509.21137].

The current literature therefore presents distributed in-memory PDHG not as a single algorithmic artifact but as a layered paradigm. At the algorithmic level, PDHG can be recast as a primal-only penalty-like method for some linear constraints, or engineered into a high-accuracy matrix-free LP solver. At the systems level, it can be mapped either to 2D distributed GPU memory with reduction-only communication or to static RRAM crossbar arrays with fixed encoded operators and exactly two in-memory MVMs per iteration. Across these variants, the recurring design principle is the same: preserve the sparse linear operator as the central computational object, keep it resident across iterations, and restructure synchronization, storage, and adaptation so that matrix-vector products remain the dominant cost.

Source: https://www.emergentmind.com/topics/distributed-in-memory-primal-dual-hybrid-gradient-pdhg