Distributed In-Memory PDHG Techniques

• Distributed in-memory PDHG is a class of first-order optimization methods that utilize distributed memory and proximal operations to solve large-scale optimization problems.
• It reformulates saddle-point problems into primal-only iterations and leverages matrix-free implementations to efficiently tackle linear programming challenges.
• The method employs adaptive step-size control and decentralized consensus strategies, enabling scalable multi-GPU and in-memory analog computing deployments.

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 (Malitsky, 2017, Applegate et al., 2021, Li et al., 12 Jan 2026, Vo et al., 25 Sep 2025).

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 (Malitsky, 2017).

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 (Malitsky, 2017).

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

$A\in \mathbb{R}^{m\times n}$0

and the saddle function

$A\in \mathbb{R}^{m\times n}$1

In this LP specialization, the dominant kernels remain sparse matrix-vector products $A\in \mathbb{R}^{m\times n}$2 and $A\in \mathbb{R}^{m\times n}$3, together with projections onto box sets, which is precisely the computational anatomy that motivates distributed-memory and in-memory deployments (Li et al., 12 Jan 2026).

2. Exact primal-only reduction and decentralized consensus

For the saddle problem

$A\in \mathbb{R}^{m\times n}$4

the specialized PDHG iteration is

$A\in \mathbb{R}^{m\times n}$5

Under $A\in \mathbb{R}^{m\times n}$6 and $A\in \mathbb{R}^{m\times n}$7, eliminating the dual recursion yields the exact primal-only form

$A\in \mathbb{R}^{m\times n}$8

and, using linearity of $A\in \mathbb{R}^{m\times n}$9,

$b\in \mathbb{R}^m$0

The sequence $b\in \mathbb{R}^m$1 generated by the primal-only method is identical to that produced by PDHG, under the simplifying initialization $b\in \mathbb{R}^m$2. The method is interpreted as a modified Tseng accelerated proximal-gradient scheme with a penalty parameter increasing to infinity; with

$b\in \mathbb{R}^m$3

it can be viewed as performing one proximal-gradient-like step for each increasing-penalty problem $b\in \mathbb{R}^m$4 (Malitsky, 2017).

The convergence theory is broader than standard saddle-point analyses. Under $b\in \mathbb{R}^m$5, $b\in \mathbb{R}^m$6, and $b\in \mathbb{R}^m$7, the averaged sequence satisfies

$b\in \mathbb{R}^m$8

more precisely

$b\in \mathbb{R}^m$9

If strong duality holds, then both $g:\mathbb{R}^n\to(-\infty,+\infty]$0 and $g:\mathbb{R}^n\to(-\infty,+\infty]$1 converge to a solution of the generalized problem, with

$g:\mathbb{R}^n\to(-\infty,+\infty]$2

Without strong duality, but with bounded solution set and $g:\mathbb{R}^n\to(-\infty,+\infty]$3 bounded below, all cluster points of $g:\mathbb{R}^n\to(-\infty,+\infty]$4 belong to $g:\mathbb{R}^n\to(-\infty,+\infty]$5, and

$g:\mathbb{R}^n\to(-\infty,+\infty]$6

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 (Malitsky, 2017).

The decentralized specialization is particularly important for distributed in-memory settings. For consensus over a connected undirected graph $g:\mathbb{R}^n\to(-\infty,+\infty]$7,

$g:\mathbb{R}^n\to(-\infty,+\infty]$8

stacking $g:\mathbb{R}^n\to(-\infty,+\infty]$9 gives

$y\in\mathbb{R}^m$0

where $y\in\mathbb{R}^m$1 is the graph Laplacian. Replacing the equivalent constraint $y\in\mathbb{R}^m$2 with $y\in\mathbb{R}^m$3 and then applying the primal-only reduction yields

$y\in\mathbb{R}^m$4

Per node,

$y\in\mathbb{R}^m$5

with a distributed implementation:

1. form $y\in\mathbb{R}^m$6;
2. exchange $y\in\mathbb{R}^m$7 with neighbors;
3. compute $y\in\mathbb{R}^m$8;
4. update $y\in\mathbb{R}^m$9;
5. update $$\min_{x\in \mathbb{R}^n}\max_{y\in\mathbb{R}^m}\quad g(x)+\langle Ax,y\rangle-\langle b,y\rangle.$$0.

Only one neighborhood communication is needed. By contrast, direct PDHG on $$\min_{x\in \mathbb{R}^n}\max_{y\in\mathbb{R}^m}\quad g(x)+\langle Ax,y\rangle-\langle b,y\rangle.$$1 requires $$\min_{x\in \mathbb{R}^n}\max_{y\in\mathbb{R}^m}\quad g(x)+\langle Ax,y\rangle-\langle b,y\rangle.$$2 and $$\min_{x\in \mathbb{R}^n}\max_{y\in\mathbb{R}^m}\quad g(x)+\langle Ax,y\rangle-\langle b,y\rangle.$$3, hence two communication rounds per iteration, and storage of $$\min_{x\in \mathbb{R}^n}\max_{y\in\mathbb{R}^m}\quad g(x)+\langle Ax,y\rangle-\langle b,y\rangle.$$4, $$\min_{x\in \mathbb{R}^n}\max_{y\in\mathbb{R}^m}\quad g(x)+\langle Ax,y\rangle-\langle b,y\rangle.$$5, and $$\min_{x\in \mathbb{R}^n}\max_{y\in\mathbb{R}^m}\quad g(x)+\langle Ax,y\rangle-\langle b,y\rangle.$$6. The primal-only scheme stores only $$\min_{x\in \mathbb{R}^n}\max_{y\in\mathbb{R}^m}\quad g(x)+\langle Ax,y\rangle-\langle b,y\rangle.$$7 and $$\min_{x\in \mathbb{R}^n}\max_{y\in\mathbb{R}^m}\quad g(x)+\langle Ax,y\rangle-\langle b,y\rangle.$$8, and the step-size condition improves from $$\min_{x\in \mathbb{R}^n}\max_{y\in\mathbb{R}^m}\quad g(x)+\langle Ax,y\rangle-\langle b,y\rangle.$$9 to $f(x)=\frac12\|Ax-b\|^2,\qquad f_*:=\min_x f(x),$0 (Malitsky, 2017).

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

$f(x)=\frac12\|Ax-b\|^2,\qquad f_*:=\min_x f(x),$1

with the associated saddle problem

$f(x)=\frac12\|Ax-b\|^2,\qquad f_*:=\min_x f(x),$2

where

$f(x)=\frac12\|Ax-b\|^2,\qquad f_*:=\min_x f(x),$3

The baseline PDHG iteration is

$f(x)=\frac12\|Ax-b\|^2,\qquad f_*:=\min_x f(x),$4

with the reparameterization

$f(x)=\frac12\|Ax-b\|^2,\qquad f_*:=\min_x f(x),$5

so the standard convergence condition becomes $f(x)=\frac12\|Ax-b\|^2,\qquad f_*:=\min_x f(x),$6 (Applegate et al., 2021).

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 $f(x)=\frac12\|Ax-b\|^2,\qquad f_*:=\min_x f(x),$7, one $f(x)=\frac12\|Ax-b\|^2,\qquad f_*:=\min_x f(x),$8 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 (Applegate et al., 2021).

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 (Applegate et al., 2021).

The convergence theory is stratified. For the LP saddle problem, baseline PDHG converges to an optimal solution if

$f(x)=\frac12\|Ax-b\|^2,\qquad f_*:=\min_x f(x),$9

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 (Applegate et al., 2021).

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 $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad x\in \arg\min f.$$0 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 $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad x\in \arg\min f.$$1, 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 $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad x\in \arg\min f.$$2 nodes, 35s at $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad x\in \arg\min f.$$3 nodes, 11 min at $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad x\in \arg\min f.$$4 nodes, and 5.4 hr at $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad x\in \arg\min f.$$5 nodes; at $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad x\in \arg\min f.$$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 (Applegate et al., 2021).

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 $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad x\in \arg\min f.$$7, 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}\quad g(x)\quad \text{s.t.}\quad x\in \arg\min f.$$8

rewritten as the saddle-point problem

$$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad x\in \arg\min f.$$9

For the Vũ–Condat form of PDHG,

$$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad A^\top A x=A^\top b.$$0

with

$$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad A^\top A x=A^\top b.$$1

the paper shows how to estimate the spectral radius of the linearized iteration from residual norms alone and to tune $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad A^\top A x=A^\top b.$$2 in order to reduce that spectral radius (Fercoq, 2024).

The central estimator is based on

$$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad A^\top A x=A^\top b.$$3

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 $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad A^\top A x=A^\top b.$$4, and local checks on the scalar sequence $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad A^\top A x=A^\top b.$$5. 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 (Fercoq, 2024).

Residual balancing remains relevant, especially as a warm-up phase. Goldstein’s rule compares primal and dual residual magnitudes and shifts the ratio $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad A^\top A x=A^\top b.$$6 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 (Fercoq, 2024).

The same paper also imposes an important limitation. For linear programs, asymptotic PDHG rate depends only on the product $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad A^\top A x=A^\top b.$$7, 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 (Fercoq, 2024).

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 (Li et al., 12 Jan 2026).

The central system idea is a 2D process or device grid with row axis $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad A^\top A x=A^\top b.$$8 and column axis $$\min_{x\in\mathbb{R}^n}\quad g(x)\quad \text{s.t.}\quad A^\top A x=A^\top b.$$9. The sparse matrix $$\min_{x \in \mathcal{X}} c^\top x \quad \text{subject to } Ax \in \mathcal{S},$$0 is partitioned into blocks $$\min_{x \in \mathcal{X}} c^\top x \quad \text{subject to } Ax \in \mathcal{S},$$1. Vectors in $$\min_{x \in \mathcal{X}} c^\top x \quad \text{subject to } Ax \in \mathcal{S},$$2, including $$\min_{x \in \mathcal{X}} c^\top x \quad \text{subject to } Ax \in \mathcal{S},$$3, $$\min_{x \in \mathcal{X}} c^\top x \quad \text{subject to } Ax \in \mathcal{S},$$4, and variable bounds, are partitioned by the column axis and replicated vertically across all devices in a given column. Vectors in $$\min_{x \in \mathcal{X}} c^\top x \quad \text{subject to } Ax \in \mathcal{S},$$5, including $$\min_{x \in \mathcal{X}} c^\top x \quad \text{subject to } Ax \in \mathcal{S},$$6 and row bounds, are partitioned by the row axis and replicated horizontally across all devices in a given row. On device $$\min_{x \in \mathcal{X}} c^\top x \quad \text{subject to } Ax \in \mathcal{S},$$7, local memory contains $$\min_{x \in \mathcal{X}} c^\top x \quad \text{subject to } Ax \in \mathcal{S},$$8, the local primal block $$\min_{x \in \mathcal{X}} c^\top x \quad \text{subject to } Ax \in \mathcal{S},$$9, $A\in \mathbb{R}^{m\times n}$00, bounds for that column block, and the local dual block $A\in \mathbb{R}^{m\times n}$01 with row bounds. This lets each device compute $A\in \mathbb{R}^{m\times n}$02, $A\in \mathbb{R}^{m\times n}$03, local primal projection, and local dual projection without gathers (Li et al., 12 Jan 2026).

One PDHG iteration has two synchronized vector reductions. First, each device computes

$A\in \mathbb{R}^{m\times n}$04

then an AllReduce over $A\in \mathbb{R}^{m\times n}$05 yields

$A\in \mathbb{R}^{m\times n}$06

followed by the local primal update

$A\in \mathbb{R}^{m\times n}$07

Second, each device forms $A\in \mathbb{R}^{m\times n}$08, computes

$A\in \mathbb{R}^{m\times n}$09

and an AllReduce over $A\in \mathbb{R}^{m\times n}$10 yields

$A\in \mathbb{R}^{m\times n}$11

followed by

$A\in \mathbb{R}^{m\times n}$12

The reflected Halpern step is fully local because current iterate and anchor are partitioned identically (Li et al., 12 Jan 2026).

Operation	Collective	Frequency
Primal step $A\in \mathbb{R}^{m\times n}$13	AllReduce over $A\in \mathbb{R}^{m\times n}$14	1/iter
Dual step $A\in \mathbb{R}^{m\times n}$15	AllReduce over $A\in \mathbb{R}^{m\times n}$16	1/iter
Halpern update	none	none
KKT evaluation	vector sums and scalar reductions	every $A\in \mathbb{R}^{m\times n}$17 iterations

This design reduces each communication to a reduction over one axis only, with communicated vector sizes approximately $A\in \mathbb{R}^{m\times n}$18 for $A\in \mathbb{R}^{m\times n}$19 and $A\in \mathbb{R}^{m\times n}$20 for $A\in \mathbb{R}^{m\times n}$21. 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 (Li et al., 12 Jan 2026).

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; $A\in \mathbb{R}^{m\times n}$22 improved from 377s to 61s; and $A\in \mathbb{R}^{m\times n}$23 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 (Li et al., 12 Jan 2026).

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,

$A\in \mathbb{R}^{m\times n}$24

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 (Vo et al., 25 Sep 2025).

The key hardware mapping is the symmetric block matrix

$A\in \mathbb{R}^{m\times n}$25

which satisfies

$A\in \mathbb{R}^{m\times n}$26

This permits one static encoded operator to serve both forward and transpose multiplies, avoiding separate programming of $A\in \mathbb{R}^{m\times n}$27 and $A\in \mathbb{R}^{m\times n}$28. 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 (Vo et al., 25 Sep 2025).

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 $A\in \mathbb{R}^{m\times n}$29, and initialization uses

$A\in \mathbb{R}^{m\times n}$30

Each iteration then uses a deterministic Nesterov-style factor

$A\in \mathbb{R}^{m\times n}$31

followed by exactly two accelerator MVMs per PDHG iteration: one for $A\in \mathbb{R}^{m\times n}$32 and one for $A\in \mathbb{R}^{m\times n}$33. The design explicitly rejects line-search-based adaptation because multiple trial MVMs per outer iteration are undesirable on IMC hardware (Vo et al., 25 Sep 2025).

Theoretical support is provided for both norm estimation and inexact PDHG under analog perturbations. For

$A\in \mathbb{R}^{m\times n}$34

the dominant eigenvalue of $A\in \mathbb{R}^{m\times n}$35 equals the dominant singular value of $A\in \mathbb{R}^{m\times n}$36,

$A\in \mathbb{R}^{m\times n}$37

which justifies Lanczos on the symmetric block matrix. Under noisy MVMs, the paper gives

$A\in \mathbb{R}^{m\times n}$38

and, for the ergodic average of the noisy PDHG iterates,

$A\in \mathbb{R}^{m\times n}$39

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 (Vo et al., 25 Sep 2025).

The RRAM evaluation uses MELISO+ as the full-stack simulation platform, with EpiRAM and TaO$A\in \mathbb{R}^{m\times n}$40-HfO$A\in \mathbb{R}^{m\times n}$41 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$A\in \mathbb{R}^{m\times n}$42-HfO$A\in \mathbb{R}^{m\times n}$43; 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$A\in \mathbb{R}^{m\times n}$44-HfO$A\in \mathbb{R}^{m\times n}$45 on the same benchmark. Accuracy is described as comparable, but not uniformly at the $A\in \mathbb{R}^{m\times n}$46 level: for gen-ip021, GPU achieved $A\in \mathbb{R}^{m\times n}$47 while both RRAM technologies were approximately $A\in \mathbb{R}^{m\times n}$48, and for gen-ip054 the RRAM gaps were $A\in \mathbb{R}^{m\times n}$49 and $A\in \mathbb{R}^{m\times n}$50. 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 (Vo et al., 25 Sep 2025).

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.