Solving linear programming in one step with RRAM-based analog matrix computing

1 Introduction

Linear programming (LP), as one of the most fundamental branches of optimization theory, plays an indispensable role in diverse applications ranging from logistics routing and power grid scheduling (Hörsch et al., 2018; Mohagheghi et al., 2018) to resource allocation and parameter optimization in distributed machine learning systems (Cai et al., 2018; Liu et al., 2018). With the rise of artificial intelligence (AI) and cloud computing, the dimensionality of LP problems in these application scenarios continues to increase (Liu et al., 2018), placing growing demands on the solving speed and energy efficiency of computational hardware.

Conventional digital solvers typically employ iterative algorithms to obtain the optimal solution of LP problems, such as the simplex method and interior-point methods (Cai et al., 2018; Vanderbei, 2020). However, these methods suffer from severe scalability issues as the problem dimensionality increases. The fundamental reason is that they inherently rely on complex matrix operations, including matrix inversion, matrix-vector multiplication (MVM), and matrix-matrix multiplication, whose computational complexity is typically polynomial, generally in the range of O n 2 to O n 3 (Wright, 1997). Under the traditional digital computing architecture, frequent data movement further limits their throughput in large-scale and real-time processing scenarios.

Analog matrix computing (AMC) technologies based on resistive random-access memory (RRAM) offer a new physical foundation for breaking through these bottlenecks (Ielmini and Pedretti, 2025). Many non-volatile memories exhibit the advantages of fast read/write speed, low read/write energy, and stable storage of single-bit or multi-bit information. Compared with other non-volatile memory devices, RRAM devices offer unique advantages, including compatibility with standard CMOS processes, simple device structure, and low fabrication cost (Ielmini and Wong, 2018; Mannocci et al., 2023; Pedretti and Ielmini, 2021). In-memory computing architectures based on multi-bit RRAM arrays eliminate the need for frequent data movement during computation acceleration (Sun et al., 2023; Wan et al., 2022). By mapping matrix entries to the conductance values of RRAM cells within the array and input vectors to input voltage vectors, RRAM-based AMC circuits can perform MVM in a single step (Hu et al., 2016; Xia et al., 2016; Yao et al., 2020). Furthermore, by constructing a global negative feedback loop, more complex matrix operations, such as matrix inversion and generalized matrix inversion, can also be accomplished in one step (Sun et al., 2019, 2020; Sun and Ielmini, 2022). Therefore, the highly parallel AMC circuits are well-suited for accelerating the complex matrix computations involved in LP solvers. Nevertheless, existing research on RRAM-based LP solvers has predominantly focused on small-scale problems (Di Marco et al., 2022; Liu et al., 2024; Vichik and Borrelli, 2014; Vichik et al., 2016) or relied on discrete iterative algorithms (Cai et al., 2018; Liu et al., 2018; Shang et al., 2020), failing to fully exploit the computing capability inherent in AMC circuits.

To address these limitations, this work proposes an AMC circuit grounded in the projection neural network (PNN) model that solves LP problems in one step. The proposed circuit directly maps the PNN dynamical model onto a global closed-loop system composed of RRAM crossbar arrays and supporting analog circuits. Compared with conventional recurrent neural network (RNN) solvers (Kennedy and Chua, 1988; Tank and Hopfield, 1986), the PNN does not require complex penalty parameter tuning and offers rigorous guarantees of global convergence (Jin et al., 2019; Liu and Zhao, 2015; Qin et al., 2017; Xia and Wang, 2000; 2004; Xia et al., 2002). These properties make the PNN-based analog solver particularly well suited for complex LP problems. To validate the proposed approach, we evaluate the circuit on assignment problems and benchmark instances from the LP Netlib test set. Circuit-level simulations show that stable solutions are obtained within approximately 15 μs for problems with up to 100 variables, achieving nearly two orders of magnitude speedup compared with the digital PNN baseline.

Furthermore, to overcome the precision limitations of AMC circuits, we introduce an analog-digital hybrid computing framework: the analog circuit rapidly generates a high-quality initial solution, which is then used to initialize a high-precision digital solver. This significantly reduces the number of required iterations, enabling an effective trade-off between speed and accuracy (Song et al., 2024; Zuo et al., 2025).

The remainder of this paper is organized as follows. Section 2 establishes the mathematical model linking LP to the PNN framework. Section 3 details the proposed circuit architecture, including the RRAM-based projection matrix computation unit and the analog neuron modules. Section 4 presents results on a simple LP problem, large-scale assignment problems, instances from the LP Netlib test set, and the proposed hybrid computing paradigm. Section 5 discusses the implications and future directions, and Section 6 concludes the paper.

2 Mathematical model of LP and PNN

A standard-form LP problem, as formulated in Equation 1, seeks to minimize a linear objective function subject to linear equality constraints and non-negativity constraints (Vanderbei, 2020; Wright, 1997):where x ∈ R n is the decision variable vector to be optimized, c ∈ R n is the cost coefficient vector, A ∈ R m × n is the constraint matrix with m < n and full row rank, and b ∈ R m is the constraint target vector. The projection matrix P = A T A A T − 1 A requires A to have full row rank. Rank deficiency of the constraint set A x = b corresponds either to redundant but consistent rows—which can be removed to yield a full-row-rank equivalent system with an unchanged feasible region—or to inconsistent rows (e.g., a relation 0 = b i with b i ≠ 0), which indicates infeasibility. As is standard in LP theory, we therefore assume without loss of generality that rank A = m, with any consistent instance preprocessed into a full-row-rank equivalent prior to circuit mapping. In practice, LP problems rarely appear directly in the standard form stated above. A general LP problem may involve objective function maximization, mixed inequality and equality constraints, or variables with non-zero lower bounds (Vanderbei, 2020). If the original problem seeks to maximize the objective function, i.e., Maximize c T x, the problem is equivalently reformulated as Minimize − c T x. For inequality constraints of the form A x ≤ b, a non-negative slack variable s ≥ 0 is appended to yield the equality A x + s = b (Dantzig and Thapa, 2003). The newly introduced auxiliary variables are augmented into the decision variable vector x. The cost coefficients associated with all slack variables are set to zero, ensuring that they do not influence the objective function. Finally, the standard form requires all decision variables to satisfy non-negativity constraints. If a variable x j is unrestricted in sign, it is decomposed into the difference of two non-negative variables: x j = x j + − x j − (Vanderbei, 2020). Through the sequential application of these transformations, any LP problem can be reduced to the standard form.

Minimize c T x Subject to A x = b , x ≥ 0 (1)

To solve the standard LP problem, we introduce the PNN model (Figure 1a). The selection of the projection neural network as the algorithmic backbone of the proposed solver is motivated by its fundamental advantages over both traditional digital iterative algorithms and earlier recurrent neural network formulations. Classical LP solvers, the simplex method and interior-point methods, are inherently sequential: the former traverses polytope vertices through dependent pivot operations with exponential worst-case complexity, while the latter requires costly matrix factorizations of O m n 2 per iteration over multiple steps (Nocedal and Wright, 2006; Wright, 1997). Earlier neural network approaches, notably the Hopfield network and its penalty-function-based variants, encode constraints through quadratic penalty terms whose coefficients must be carefully tuned for each problem instance (Kennedy and Chua, 1988; Tank and Hopfield, 1986). Unlike early Hopfield networks or penalty-function-based RNNs, the PNN utilizes a projection operator P Ω · to enforce the state trajectory of the system within the feasible region. For example, the projection operator P Ω x = argmin z ∈ Ω z − x identifies the projected point in the feasible region nearest to the current variable x, thereby avoiding convergence problems caused by improper selection of penalty parameters. According to PNN theory, the state evolution equation for solving the LP problem can be described as:where y ∈ R n is the neuron state variable vector, and y + = max 0 , y is the output after the rectified linear unit (ReLU) activation function, representing the current optimization solution x. The parameter ϵ is a time constant that determines the convergence speed of the system. The matrix P is the projection matrix, defined as P = A T A A T − 1 A, which projects any vector onto the row space of the constraint matrix A. Correspondingly, I − P projects vectors onto the null space of A, ensuring that gradient updates do not violate the equality constraint A x = b. The bias term u is defined as − I − P c + A T A A T − 1 b, which incorporates the projected component of the objective function gradient c and the correction component required to satisfy the constraint b (Jin et al., 2019; Liu and Zhao, 2015; Qin et al., 2017).

d y d t = 1 ϵ − I − P y + I − 2 P y + + u (2)

The equilibrium point of this differential equation corresponds to the optimal solution of the original LP problem. Existing literature has rigorously proven the global asymptotic stability of this dynamical system in the Lyapunov sense, providing a solid theoretical guarantee for the reliability of hardware implementation (Jin et al., 2019; Liu and Zhao, 2015; Qin et al., 2017; Xia and Wang, 2000; 2004; Xia et al., 2002). The correspondence between the equilibrium of Equation 2 and the optimum of the standard-form LP can be established directly from the Karush–Kuhn–Tucker (KKT) conditions. At an equilibrium, decompose y = y + − y − with y + = max ⁡ 0 , y, y − = max ⁡ 0 , − y, and use the idempotence P 2 = P. Setting d y / d t = 0 in Equation 2 gives 0 = − P y + + I − P y − − I − P c + A T A A T − 1 b . Left-multiplying by P, and using P I − P = 0 and P A T = A T, yields P y + = A T A A T − 1 b; applying A then gives primal feasibility A y + = b. Left-multiplying by I − P gives I − P c − y − = 0, i.e., c = A T λ + s with dual variable λ = A A T − 1 A c − y − and dual slack s = y − ≥ 0, which is dual feasibility. Since y i + y i − = 0 component-wise, s T x = 0, which is complementary slackness. The equilibrium therefore satisfies the full KKT system, so that x * = y + is the LP optimum and the dual solution is recoverable from the circuit.

The global asymptotic stability of this dynamical system has been rigorously established in the PNN literature (Jin et al., 2019; Liu and Zhao, 2015; Qin et al., 2017; Xia and Wang, 2000; 2004; Xia et al., 2002), under the assumption that A has full row rank so that A A T is invertible and P is well defined, and that the problem is feasible and bounded. The proof applies unchanged to the scaled circuit equation, since the optimum-preserving voltage and conductance scalings introduced in Section 3 leave the location of the equilibrium invariant. For the non-ideal circuit (finite OA gain, offset, conductance error, parasitic resistance) these guarantees hold approximately; the resulting deviation of the equilibrium is characterized empirically in Section 4 and Section 5.

3 Architecture of RRAM-based AMC LP solver

In order to translate the above matrix differential equation into a physical circuit that solves LP in one step, we propose a hardware architecture integrating a closed-loop AMC unit with analog neuron modules. As illustrated in Figure 1b, the proposed analog LP solver comprises two main components. The first is the Projection Unit, built on RRAM crossbar arrays and functional peripheral circuits, which completes the complex matrix operations associated with the projection matrix P in one step. The second is the Analog Neuron Module, composed of operational amplifiers (OAs), which realizes the integration, addition, subtraction, and nonlinear activation required by the differential equation. These two components are hardwired into a global negative feedback system. Once the circuit is powered on and voltage signals encoding vectors b and c are applied, the system settles within microseconds, yielding the optimal solution y + at the neuron output terminals—completing the LP solution in one step.

Throughout this paper, one step refers to two specific properties: (i) the composite matrix operation P − I v i n + A T A A T − 1 b, which costs O m n 2 per evaluation in the digital domain, is executed by the Projection Unit as a single physical operation whose latency is set by circuit settling rather than arithmetic complexity; and (ii) the LP solution is obtained in a single continuous-time transient after input application, without discrete iterations, conductance reprogramming, or external control flow. It does not mean an instantaneous algebraic solution: the closed-loop system converges with a time constant determined by ϵ and the loop dynamics. The proposed solver targets standard-form LP instances that are feasible, bounded, and reducible to a full-row-rank equality system; coefficient scaling into the admissible voltage and conductance windows, reduction to standard form, and elimination of consistent redundant constraints are assumed to be performed as preprocessing prior to mapping the problem onto the hardware. The architecture does not, in its present form, incorporate in-circuit detection of infeasibility or unboundedness, nor does it resolve degenerate optima, multiple optima, or failure of strict complementarity. It is accordingly presented as a specialized accelerator as a proof-of-concept for well-conditioned, feasible, bounded standard-form LPs, rather than as a general-purpose LP solver.

3.1 Projection unit

Computing the projection matrix P = A T A A T − 1 A in Equation 2 entails matrix multiplication and matrix inversion, operations that are prohibitively expensive in the digital domain. As shown in Figure 1b, this work realizes this computation through a topological extension of an RRAM-based analog matrix generalized inverse circuit (Sun and Ielmini, 2022). The unit comprises two RRAM crossbar arrays storing the transpose of the constraint matrix, namely A T. Each element A i j is mapped to the conductance value G i j of RRAM. Each RRAM cell is modeled as an ohmic resistor operating within its linear conduction regime; provided this condition holds, both 1R and 1T1R configurations admit a linear resistive equivalent. The aggregate effect of device-level non-idealities—conductance programming error, cycle-to-cycle and device-to-device variability, drift, relaxation, read noise, and read disturb—is approximated by a single Gaussian-distributed deviation of the programmed conductance relative to the maximum conductance state, the impact of which on solution accuracy is analyzed in Section 5. Since RRAM conductances are inherently non-negative, matrices containing negative elements are accommodated via a column-line differential structure or a bias conductance scheme, as shown in Supplementary Figures S2, S3 (Luo et al., 2025).

Within the Projection Unit, the column lines of the right array are respectively connected to the non-inverting inputs of a group of OAs, whereas the column lines of the left array are connected to the outputs. Exploiting the virtual short and virtual open properties of the OAs, the matrix equation realized by the right-side array is:where k is a scaling factor and v o u t is the row-line voltage vector of the right-side array, serving as the module output. A group of transimpedance amplifiers (TIAs) is incorporated, whose inputs are connected to the row lines of the left array and outputs are connected to the row lines of the right array. As a result, the left array essentially accomplishes the MVM function. The matrix equation realized by the left-side array is:where v i n is the module input voltage vector, applied to the row lines of the left-side array via unit conductances, and v r is the column-line voltage vector of the left-side array. Substituting Equation 4 into Equation 3 yields Equation 5:

b k + A · v o u t = 0 (3)

k v o u t = A T v r − v i n (4)

b k + 1 k A A T v r − A v i n = 0

A A T v r = A v i n − b

v r = A A T − 1 A v i n − b (5)

Back-substitution gives:

v o u t = 1 k A T A A T − 1 A v i n − b − v i n

v o u t = 1 k P − I v i n − 1 k A T A A T − 1 b (6)

This derivation confirms that the Projection Unit executes the full set of matrix operations required by the PNN in a single physical step. Furthermore, the Projection Unit is capable of functioning independently. When the input vector b is set to zero, the output of the Projection Unit becomes v o u t = 1 k P − I v i n. In applications beyond the LP solver, the Projection Unit can be employed to accelerate any computational task involving the projection operator.

3.2 Analog neuron module

To complete the hardware realization of the PNN differential equation, each state variable y j is served by a dedicated analog neuron circuit containing O A 1, O A 2, and O A 3, as illustrated in Figure 1b. The OAs are described by a macromodel with supply rails of ± 5 V, an input common-mode range of ± 4.5 V, an output swing of ± 5 V, a slew rate of 22 V / μ s, a − 3 dB bandwidth of 16 MHz, an input offset of ± 0.8 mV, and a settling time of 350 ns. The complete circuit parameters, including the positive-domain, real-domain, and conductance-compensated real-domain implementations, are given in Supplementary Figures S1–S3. O A 1 constitutes the core integrator of the analog neuron. By exploiting the charge–discharge characteristics of a feedback capacitor, it performs continuous-time integration, thereby realizing the ϵ d y / d t term in Equation 2, where ϵ = C f / G 0 is the time constant determined by the feedback capacitance and conductance. Its inputs originate from the Projection Unit output, the outputs of O A 2 and O A 3, and an external bias voltage encoding the cost vector − c. By Kirchhoff’s current law, the voltages at the non-inverting and inverting input nodes of the j-th OA₁ are obtained from the respective resistive networks in Equation 7 as follows:

v O A 1 , j + = − c j · G 0 + y j + · G 0 G 0 + G 0 + k + 1 G 0 = − c j + y j + k + 3

v O A 1 , j − = y j + ϵ d y j d t · G 0 + v o u t , j · k G 0 + v i n , j · G 0 G 0 + k G 0 + G 0 + G 0 = y j + ϵ d y j d t + k v o u t , j + v i n , j k + 3 (7)

Invoking the virtual-short condition v O A 1 , j + = v O A 1 , j − and selecting the resistors and capacitors to match the feedback parameter k of the Projection Unit, the common denominator k + 3 cancels. Consequently, expressing the outputs y j of all the O A 1 in vector form, the output voltage y directly represents the state variable

y + ϵ d y d t = − c + y + − k v o u t + v i n (8)

O A 2 implements the nonlinear activation y + = max 0 , y required by the PNN. A single-supply voltage follower configuration truncates negative voltages: when y < 0, the output is clamped to 0; when y > 0, the output tracks linearly. O A 3 functions as an analog subtractor whose output is routed to the input of the Projection Unit, thereby closing the global feedback loop of the solver. The voltages at the non-inverting and inverting input nodes of the j-th OA₃ are expressed in Equation 9:

v O A 3 , j + = − c j · G 0 + y j + · 2 G 0 G 0 + 2 G 0 + G 0 = − c j + 2 y j + 4

v O A 3 , j − = v i n , j · G 0 + y j · G 0 G 0 + G 0 + 2 G 0 = v i n , j + y j 4 (9)

Invoking the virtual-short condition v O A 3 , j + = v O A 3 , j −, the output voltage vector of O A 3 is given by:

v i n = − c + 2 y + − y (10)

Substituting Equation 6 of the projection module into Equation 8 yields Equation 11:

y + ϵ d y d t = − c + y + − P v i n − A T A A T − 1 b

ϵ d y d t = − P v i n − y − c + y + + A T A A T − 1 b (11)

Finally, inserting the subtractor relation of Equation 10 and collecting terms in y, y +, and c yields the governing differential equation of the analog solver:

ϵ d y d t = − I − P y + I − 2 P y + − I − P c + A T A A T − 1 b (12)

Equation 12 is identical in form to the canonical PNN dynamics, confirming that the three-amplifier neuron, in conjunction with the Projection Unit, faithfully maps the PNN differential equation onto hardware. Although the projection-unit output in Equation 6 contains the term − A T A A T − 1 b, the inverting configuration of the corresponding OA stage within the neuron restores the + A T A A T − 1 b sign at the summing node, confirming that the closed loop implements the canonical PNN rather than a sign-scaled variant. The solution is read at y + when the rectified output changes by less than a prescribed tolerance over a fixed observation window. Coefficient normalization is performed only through optimum-preserving similarity scalings: scaling a constraint row of A and the corresponding entry of b by a common positive factor leaves the feasible set unchanged, and a uniform positive scaling of c leaves the optimizer unchanged. Negative matrix entries are encoded as a differential pair of conductances, with a true zero realized by equal conductances; negative entries of b and c are applied directly as negative reference voltages.

4 Results

To validate the proposed one-step analog LP solver, a series of experiments of increasing complexity were conducted. A small-scale LP problem was first constructed as a proof-of-concept demonstration. As shown in Figure 2a, by introducing slack variables, the original problem is converted into a standard LP problem with five variables and three equality constraints. Since the constraint matrix A contains real-valued entries, it is mapped onto the conductance values of the RRAM array through a matrix-splitting scheme G A = G B − G C, where G B encodes the positive elements of A, and G C encodes the absolute values of the negative elements of A, shown in Figure 2b. The objective vector c and constraint vector b are applied as voltage-source inputs.

As shown in Figures 2c,d, transient responses from numerical simulation of the PNN-based digital solver and circuit simulation with Simulation Program with Integrated Circuit Emphasis (SPICE) are compared. The digital solver was implemented in MATLAB using the same PNN model and executed on an Intel® Core™ i5-14600KF processor. The digital solver shows convergence to the optimal solution (0.3000, 0.6000, 0.0000, 0.0000, 0.4500) after 300 iterations. Although the problem size is relatively small, the digital solver still requires approximately 300 μs to reach the final solution, owing to the computationally intensive operations associated with the projection matrix.

The proposed analog LP solver produces steady-state output voltages of approximately (0.3003, 0.6069, −0.0046, −0.0052, 0.4564) within a few microseconds, demonstrating excellent agreement. The plotted analog outputs are the rectified neuron states y + read after the ReLU stage and are identified with the optimization solution x. The small negative values observed for variables whose expected values are zero (e.g., −0.0052 V) do not reflect a failure of the rectification; they arise from O A 2 operating far from its common-mode point, where the finite open-loop gain produces a residual offset of a few millivolts about zero (Rao et al., 2023; Zhou et al., 2025). For LP problems with moderate accuracy requirements, this residual is small relative to the dynamic range of the solution variables and does not materially affect the validity of the result. The corresponding raw, unrectified states y are reported in Supplementary Figure S4.

Figures 2e,f depict the solving trajectories in the phase plane, color-coded by iteration index and by simulation time, respectively. Figure 3 traces the evolution of the objective function value and constraint residual ∥ A x − b ∥ for both the circuit and the numerical simulation, where · denotes the l 2-norm (Euclidean norm) of a vector. As shown in Figure 3a, the objective value of the digital solver rapidly stabilizes, while the constraint residual is progressively refined and eventually converges to a high precision after approximately 600 iterations. The simulation results of the analog LP solver presented in Figure 3b indicate that it could converge to a near-identical objective value. However, due to the parasitic effects and non-idealities present in the circuit, the constraint residual can only reach a relatively low precision when the circuit output stabilizes.

4.1 Assignment problems

To further demonstrate the capability of the proposed analog LP solver, the problem size was scaled to classic assignment problems involving 64 and 100 variables. Taking the assignment of tasks to workers as an example, as illustrated in Figure 4a, the assignment problem can be interpreted as assigning n workers to n jobs, where each worker must be assigned to exactly one job and each job must be assigned to exactly one worker. The element c i j in the cost matrix C represents the cost, time, or resource consumption associated with assigning worker i to job j. The goal is to determine an assignment matrix X = x i j, where x i j = 1 indicates that worker i is assigned to job j, and x i j = 0 otherwise, such that the total assignment cost is minimized (Liu and Zhao, 2015). Equivalently, the problem seeks to identify n entries in an n × n cost matrix, one from each row and one from each column, that minimize the total cost. Accordingly, an 8 × 8 matrix results in 64 decision variables, whereas a 10× 10 matrix yields 100 decision variables. The corresponding standard LP formulation is given in Equation 13:

Minimize ∑ i = 1 n ∑ j = 1 n c i j x i j Subject to ∑ j = 1 n x i j = 1 , ∑ i = 1 n x i j = 1 , x i j ≥ 0 (13)

In the assignment problem, the row- and column-selection constraints yield, after eliminating one redundant constraint (the last-column constraint, owing to the redundancy of the constraint set), a reduced full-row-rank constraint matrix of size 2 n − 1 × n 2 with an all-ones right-hand-side vector of dimension 2 n − 1; the explicit construction follows (Liu and Zhao, 2015) and is given in Supplementary Data S1. The cost matrices used for the 8 × 8 and 10× 10 instances are listed in Supplementary Tables S1, S2. Figure 4b illustrates the constraint matrix mapped onto the conductance values of the RRAM array. As shown in Figures 4c,d, the 64-variable assignment problem results are compared between digital and analog LP solvers. The RRAM-based analog LP solver reaches steady state within approximately 15 μs, whereas the digital PNN baseline requires roughly 1.9 ms, corresponding to nearly two orders of magnitude speedup relative to the digital PNN reference on this problem size. This comparison is made against the same PNN model executed digitally and is not a claim of superiority over optimized commercial or open-source LP solvers. Figure 5 illustrates the cost function vector c in matrix form for the 100-variable assignment problem. Figure 5a presents the theoretical optimal solution obtained by the digital solver, while Figure 5b shows the optimal solution produced by the analog LP solver. Owing to the increased problem size, non-idealities in the circuit cause the outputs to deviate from exact 0 and 1 values. For instance, in the element selection of the 6th column, the analog LP solver exhibits a 30% confidence in selecting the element in the 6th row and a 70% confidence in selecting the element in the 7th row. Here, “confidence” denotes the normalized analog output voltage of each variable. A 50% threshold is applied solely as a post-processing step specific to the binary assignment problem, whose polytope vertices are binary and which exhibits a sufficient margin between selected and unselected entries. It is not a general procedure: for LPs with continuous optima, thresholding can destroy both feasibility and optimality, and is therefore not applied to the continuous LP cases. By applying this threshold-based correction to the output, the analog LP solver is still able to obtain an optimal solution consistent with that of the digital PNN baseline.

Beyond solving speed, energy efficiency is a critical metric for practical deployment of analog optimization hardware. For an LP problem with an m × n constraint matrix A, the proposed LP solver requires n TIAs and m OAs in the Projection Unit, and 3 n OAs in the Analog Neurons. If the real-domain RRAM matrix is implemented using analog inverters, additional n + m OAs are required, whereas the real-domain RRAM matrix based on the conductance-compensation method introduces no additional OAs. In addition, n + m DACs are required to provide the input voltages, and n ADCs are needed to read out the results. The power consumption of the OAs and ADCs is on the order of hundreds of μW, that of the DACs is on the order of μW, while the power consumed by the RRAM array is nearly negligible (Ielmini and Wong, 2018; Mannocci et al., 2023; Pedretti and Ielmini, 2021; Sharma and Nath, 2024; Zuo et al., 2025). Consequently, for a 100-variable assignment problem, the total circuit power consumption is approximately 55.94 mW (67.84 mW considering real-domain RRAM matrix using analog inverters). Given that the circuit reaches steady state within 15 μs, the corresponding computation energy is only 0.8391 μJ, according to our previous estimation method (Zuo et al., 2025). For writing the RRAM array, programming a single cell requires approximately 0.4 μW. Under a parallel-write configuration, the write-verify method enables an entire row or column of RRAM cells to be simultaneously programmed to their precise target conductance states within a few tens of pulses. Consequently, when the constraint matrix of the LP problem needs to be refreshed, approximately 38 μs is required to write the matrix, corresponding to an energy of 0.1155 μJ. These power and energy estimates account for the OAs, DACs, and ADCs; the contributions of reference generation, control logic, and data-handling blocks lie outside the present analysis, so these are not a full system-level accounting. The estimate further refers to the per-solve computation energy; the one-time energy for programming the constraint matrix A onto the RRAM array for each new instance is accounted for separately. The reported speedup and energy comparison are made against the digital PNN baseline. To estimate the energy consumption of the digital baseline, we measured the runtime power of an Intel i5-14600KF CPU using the Intel Extreme Tuning Utility (Intel XTU) while running the LP solver implemented in MATLAB. The CPU consumed approximately 24.63 W during execution; combined with a runtime of 1.9 ms, the total energy dissipation reached approximately 46.8 mJ, which is significantly higher than that of the proposed analog LP solver (Li et al., 2025).

4.2 Standard test set evaluation

The solver is further benchmarked on the SC50A and SC50B instances from the LP Netlib test collection. Both problems belong to the SC (staircase-structured) family of benchmark instances (Koch, 2004). Staircase-structured LPs commonly arise in multi-period production planning, inventory management, and dynamic resource allocation scenarios, where decisions in one time period are coupled to adjacent periods through linking constraints, giving the constraint matrix a characteristic block-staircase sparsity pattern. As shown in Figure 6a, SC50A comprises 50 rows and 78 columns with 160 non-zeros, while SC50B has the same row–column dimensions but with 148 non-zeros (Figure 6b). The dimensions reported here refer to the transformed standard form rather than the original Netlib listing: the original SC50A and SC50B instances comprise 50 inequality/equality constraints over 48 variables, which become a 50-row, 78-variable equality system after appending slack variables to convert the inequality constraints to the standard form A x = b, x ≥ 0. The constraint vector b is scaled down by a factor of 200 to map the problem into the admissible voltage range. The exact transformed A, b, and c are provided in Supplementary Data S2. Notably, the cost coefficients in both SC50A and SC50B were negated when included in the Netlib repository, meaning the original problems are maximization problems that have been converted to the standard minimization form. Despite their modest size, these instances are widely used as standard validation benchmarks for LP solvers because their staircase coupling structure and near-degenerate vertex configurations can expose numerical sensitivity and convergence weaknesses in both simplex and interior-point methods (Gleixner et al., 2016).

The solution accuracy of the analog LP solver is evaluated by comparing the circuit output with the digital numerical output for every component of the solution vector. For SC50A (Figure 6c), the data points cluster tightly along the diagonal, with only minor deviations observed for a few slack variables (c i = 0) in the intermediate range. The decision variable (c i = − 1) falls almost exactly on the diagonal, indicating that the circuit correctly identifies the optimal value of the objective-relevant component. For SC50B (Figure 6d), an even tighter clustering along the diagonal is observed, demonstrating excellent agreement with the numerical reference solution. For SC50A and SC50B, respectively, the analog solver attains objective values c T x = −0.3109 and −0.3457 (digital reference: -0.3229 and −0.3500) with primal residuals ∥ A x − b ∥ = 0.0215 and 0.0279; the full solution vectors are listed in Supplementary Tables S3, S4. In both cases the output error relative to the digital solution remains acceptably low, confirming that the analog LP solver reliably handles real-world-origin benchmark instances with staircase coupling structure while preserving one-step analog computing.

4.3 Analog-digital hybrid computing paradigm

While the one-step AMC solver delivers substantial speed advantages, certain applications demand higher numerical precision than analog circuits alone can guarantee (Zuo et al., 2025). To reconcile analog speed with digital accuracy, we propose a hybrid computing paradigm. The analog LP solver first produces a near-optimal approximate solution on a microsecond time scale. This solution then serves as a warm-start seed for a high-precision digital iterative solver. Conventional digital iterative algorithms are highly sensitive to initialization quality. A distant initial point forces the algorithm through a protracted search before it enters the convergence basin. The seed furnished by the one-step analog LP solver, by contrast, already resides deep within the attraction basin of the global optimum. The digital solver used in the hybrid experiment is the same PNN discretization employed throughout this work (step size 0.1, zero initialization, terminated at a prescribed maximum iteration count); the digital refinement repairs any residual infeasibility of the analog seed, which is not assumed to satisfy A x = b or x ≥ 0 exactly. As shown in Figure 7, seeding the digital solver with the analog output reduces the iteration count by up to 54.4% relative to random initialization for equivalent constraint satisfaction. The iteration-reduction behavior across 100 instances is shown in Supplementary Figure S6. To reach constraint-residual thresholds of 10 − 4, 10 − 6, and 10 − 8, warm-starting reduces the mean iteration count by 15.4%–37.9%. Because the analog computation time is negligible compared with even a single digital iteration, total system solving time is substantially reduced. This hybrid strategy harnesses rapid-convergence strengths of analog matrix computing to offset the slow initial-phase behavior of digital methods, charting a compelling direction for next-generation high-performance optimization architectures.

5 Discussion

This work has presented a systematic study of one-step LP solving through RRAM-based analog matrix computing, spanning mathematical modeling, circuit architecture, and benchmarking. The principal advantage of the proposed solver is that, owing to the intrinsic parallelism of RRAM crossbar arrays, the analog LP solver based on the PNN model can solve LP problems in one step. However, several aspects warrant further discussion.

Regarding precision, the analog solver inevitably introduces errors from device-level non-idealities, RRAM conductance variation and parasitic wire resistance. To quantify parasitic wire resistance, the 64-variable assignment problem was re-simulated with parasitic line resistances of 0.5, 1, and 2 Ω (Supplementary Figure S7). A slight IR-drop emerges at the output of the analog LP solver once the wire resistance reaches 2 Ω. The influence of wire resistance is therefore expected to become more pronounced when larger-scale and denser matrices are mapped onto the RRAM array. This effect can be mitigated by optimizing the array interconnect layout to shorten and widen lines, operating the RRAM devices in a higher-resistance state to increase the device-to-parasitic ratio, and refining the write strategy to compensate for systematic conductance errors. Moreover, RRAM devices are subject to stochastic write variability, read noise, and conductance drift (Ielmini and Wong, 2018; Mannocci et al., 2023; Pedretti and Ielmini, 2021; Rao et al., 2023; Zhou et al., 2025), all of which perturb the physically encoded projection matrix and may degrade solution quality. It is therefore essential to evaluate how gracefully the solver’s performance degrades as device-level non-idealities intensify. To quantify this effect, the 64-variable assignment problem was re-simulated with RRAM conductance variation. The device conductance variation is modeled using a Gaussian distribution (Pedretti et al., 2021), whose mean corresponds to the target conductance state and standard deviation varies from 0% to 6% of the maximum conductance value. As shown in Figure 8, the analog LP solver maintains correct outputs, recovering the optimal assignment after thresholding, for conductance variations up to 5%. Beyond this level, accumulated conductance errors propagate through the Projection Unit and begin to produce erroneous assignment decisions. These results indicate that the proposed architecture possesses appreciable robustness against realistic device-level variability, and that the binary nature of assignment solutions provides a natural error-tolerance margin through the thresholding operation. We note that, for continuous LP problems such as SC50A and SC50B, conductance deviation affects solution quality substantially more than for the assignment problem, where thresholding masks small analog errors. The accuracy requirement on the devices is consequently demanding relative to most present-day RRAM; with continued improvement in device precision (Rao et al., 2023; Song et al., 2024), one-step analog solution of continuous LP under these non-idealities is expected to become increasingly attainable.

It is also instructive to position the proposed solver against prior analog and memristor-based optimization circuits, as shown in Table 1. Early analog LP/QP circuits, such as the Tank-Hopfield network (Tank and Hopfield, 1986), establish that optimization dynamics can be embedded directly in continuous-time hardware, but they rely on penalty-function formulations whose parameters must be tuned per instance and were demonstrated only on small problems. Although the analog LP/QP circuit proposed by Vichik and Borrelli (Vichik and Borrelli, 2014) was demonstrated on relatively large-scale problems, accurate solutions require an a priori tuning of the cost function U cost, which must be determined through trial and error. Liu carried out silicon verification of the analog LP/QP circuit proposed by Vichik; however, it was tested only on relatively small-scale problems (Liu et al., 2024). The memristor neural networks for LP/QP of Di Marco et al. (2022) introduce resistive devices into the optimization loop, yet the dynamical models mapped onto the hardware remain comparatively simple and the validated problem sizes are small. The interior-point accelerator of Cai et al. (2018) instead employs the memristor crossbar only to accelerate the matrix-vector multiplications within each iteration of a digital algorithm, so that solving a single LP instance still requires many sequential iterations and external digital control. Shang et al. (2020) refined the iterative procedure for LP solving by means of a Douglas-Rachford framework; nevertheless, for practical instances the method may still be limited by the large iteration count and the latency associated with data movement. In contrast, the present work maps the full PNN dynamical system onto a closed-loop RRAM Projection Unit and analog neurons, so that the complete LP solution, including the costly projection-matrix computation A T A A T − 1 A, which is O m n 2 in the digital domain, is obtained in a single continuous-time transient without discrete iterations. With respect to speed and scale, our solver reaches stable solutions within approximately 15 μs for problems with up to 100 variables and is validated not only on synthetic instances but also on the real-world-origin Netlib SC50A and SC50B benchmarks, a problem scale and benchmark scope wider than those reported in the aforementioned analog LP studies. With respect to power, the full peripheral interface (OAs, DACs, and ADCs) consumes approximately 55.94 mW, corresponding to a per-solve computation energy of approximately 0.8391 μJ. With respect to accuracy, the analog solver attains moderate precision (a primal residual on the order of 10 − 2) directly, which is sufficient for binary problems such as assignment after thresholding and which can be refined to high precision through the analog-seed/digital-iteration hybrid paradigm. We therefore regard the principal advance of this work as the architecture-level mapping of a PNN LP dynamical system onto analog matrix-computing hardware, together with its accuracy and robustness characterization, rather than as a claim of novelty over analog matrix inversion or memristor-assisted optimization in general.

TABLE 1. Comparison of the proposed solver with prior analog and memristor-based LP circuits.

Finally, the proposed framework generalizes naturally to broader classes of convex optimization problems, including quadratic programming, which is central to model predictive control and support vector machine training (Liu et al., 2024). The PNN formalism accommodates these problem classes with only minor modifications to the RRAM-based projection matrix computation unit, suggesting a versatile analog matrix computing platform for one-step optimization across diverse domains.

6 Conclusion

We have proposed and validated an RRAM-based AMC circuit that solves LP problems in one step by leveraging the PNN model. The architecture couples a Projection Unit built on RRAM crossbar arrays and analog neuron modules to form a closed-loop dynamical system that converges to optimal solutions within microseconds. Simulations on problems ranging from 5 to 100 variables and on the Netlib SC50A and SC50B benchmarks demonstrate approximately two orders of magnitude speedup over the digital PNN baseline, with accuracy sufficient for practical applications. The accompanying analog-seed/digital-iteration hybrid paradigm bridges the gap between analog speed and digital precision, cutting digital iteration counts by up to 54.4%. These results establish one-step RRAM-based AMC as a promising platform for real-time, energy-efficient optimization at scale.

Correction note

This article has been corrected with minor changes. These changes do not impact the scientific content of the article.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnano.2026.1882288/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Sajid Husain, University of California, Berkeley, United States

Reviewed by: Muhammd Umair Khan, Khalifa University, United Arab Emirates

Yang Li, Peng Cheng Laboratory, China