Memristive-Friendly Hadamard Reservoir Computing: Structured, Multiplier-Free Recurrences at Scale

Motivation

• The structured recurrence addresses dense ESN costs of O(N^2) weights, operations, and physical conductances while targeting reservoir-size regimes where predictive quality improves with N.Minimal cycle reservoirs are cheap but propagate information one position per step, whereas dense orthogonal recurrences improve memory at dense-matrix cost.Source: §1 · pp. 2–3

Model contribution

• Hadamard reservoirs replace dense coupling with M = D2 bHN Π D1, using O(N) parameters and O(N log N) additions and subtractions without materializing a matrix.The operator combines sign diagonals, a permutation, and a normalized Walsh-Hadamard transform; a single binary input connection per unit removes analog-valued input matrices.Source: §1 · pp. 3–4
• The H-ESN and MF-H-ESN retain standard or memristive neuron dynamics while replacing only recurrent coupling, combining exact orthogonality, global one-step mixing, and multiplier-free implementation.The MF-H-ESN reduces recurrence to sign flips, fixed routing, and pairwise additions while preserving memristive-device dynamics.Source: §1 · p. 3

Theory

• The analysis derives an echo state condition tight in recurrent scaling ρM and a design-time-fixed noise gain, while separating the roles of randomization elements through spectral and mixing analysis.Exact orthogonality provides uniform memory and predictable responses to runtime perturbations.Source: §1 · pp. 3–4

Hardware evaluation

• The study measures recurrent-step wall-clock time on three hardware platforms and memory, showing device-dependent dense-product crossover but hardware-independent memory advantage.The paper evaluates computational and storage costs as distinct consequences of the structured operator.Source: §1 · p. 4

Empirical evaluation

• The models are evaluated on twenty classification and seven regression benchmarks against dense, cycle, dense-orthogonal, and fully trainable gated recurrent baselines up to N = 8192.The evaluation includes critical-difference analysis.Source: §1 · p. 4

Ablations and robustness

• Ablations show that permutation protects memory near the stability edge, whereas sign diagonals provide representational richness in nonlinear states.The robustness study measures accuracy costs from state noise, in-loop quantization, device mismatch, and faults in signs and routing.Source: §1 · p. 4

2.1. Echo State Networks

• ESNs use an untrained recurrent nonlinear layer with a trainable linear readout, while the ESP makes reservoir states depend on input history rather than initial conditions.The leaky formulation uses recurrent and input weights, bias, tanh nonlinearity, and leak rate α; spectral radius ρ is commonly tuned below one, although ||Wh||2 < 1 is sufficient for every input.Source: §2.1 · p. 5

2.2. Memristive-friendly Reservoir Computing

• The memristive-friendly model derives neuron dynamics from potentiation-depression kinetics, with voltage-window rescaling and state updates governed by device conductance changes.The retention factor γ −εS(zt) makes volatility γ and integration step ε jointly analogous to the ESN leak rate α, while s controls rescaling steepness.Source: §2.2 · p. 6

2.3. Structured and hardware-friendly reservoir topologies

• Hardware-oriented reservoirs range from deterministic cycle topologies to dense orthogonal matrices [13, 14], but structured sign-permutation-Walsh-Hadamard operators offer a middle ground at O(N log N) cost.Dense orthogonal matrices maximize short-term memory but retain dense storage, application, and O(N^3) factorization costs; the proposed operator is used recurrently in standard and memristive-friendly ESNs, emphasizing hardware realizability, stability, and scaling.Source: §2.3 · p. 7
• Unlike cycle reservoirs’ periodicity and Fourier-basis characterization, the proposed operator provides global mixing in a single recurrence application.This comparison motivates analyzing the consequences of replacing cyclic structure with global mixing in Section 5.Source: §2.3 · p. 7

3.1. Structured orthogonal recurrence

• The proposed operator M = D2 bHNΠD1 is exactly orthogonal, mixes the full state, and uses a single scalar ρM to control recurrent strength.Its factors provide sign randomization, permutation, and Walsh-Hadamard mixing while preserving norm and giving every eigenvalue of ρMM modulus ρM.Source: §3.1 · p. 8, §3.2 · p. 9

3.3. Cost analysis

• The recurrent step costs O(N log N) additions and subtractions and stores O(N) discrete parameters instead of O(N^2) dense weights, without materializing the matrix.The representation uses 2N signs, N routing indices, and one scaling constant; its total description requires O(N log N) bits.Source: Figure 1 · p. 9, §3.3 · p. 10

3.2. Model definitions

• The H-ESN and MF-H-ESN replace dense recurrence with the same structured operator, while the latter preserves the original potentiation-depression kinetics and leaky integration.The operator changes only how the previous state is mixed before the memristive nonlinearity.Source: §3.2 · p. 9, Figure 2 · p. 10
• Both reservoirs also support a simplified input variant with one randomly selected binary-weighted input connection per unit, reducing input wiring from Nd to N physical connections.The simplified pathway removes analog-valued input and bias parameters through deterministic binary scalings and a global gain.Source: §3.2 · p. 10, §3.3 · p. 11

3.4. Hardware considerations

• The structured recurrence enables memristive hardware based on interconnects, polarity-controlled summations, and simple scalings, while circuit-level butterfly realization remains future work.The design retains global mixing but this section studies algorithmic consequences rather than presenting a physical circuit.Source: §3.4 · p. 11

4.1. Wall-clock time across hardware platforms

• The structured recurrent step requires O(N log N) additions and subtractions instead of O(N^2) dense multiply-accumulates.Timing uses fp32 batches of 256 states across an Apple M3 Pro CPU, NVIDIA T4 GPU, and NVIDIA A100 GPU.Source: §4.1 · pp. 11–12
• At N = 32768, the structured recurrence reaches about 50× speedup on the Apple M3 Pro CPU and NVIDIA T4, while the A100 reaches 3.5× at N = 16384.The crossover occurs near N ≈1024 on the CPU and T4 but near N ≈8192 on the A100, where the structured step remains within about ten percent below that size.Source: §4.1 · p. 12, Figure 3 · p. 13

4.2. Memory footprint

• At N = 8192, structured storage is 15.4 kB versus 268 MB for dense coupling, a reduction reaching 1.7 × 10^4.Dense coupling stores 32N^2 bits and 6.7 × 10^7 analog values, whereas the structured operator uses 2N signs, N routing indices, and one gain.Source: §4.2 · p. 14
• Structured coupling reduces programmable analog values from N^2 to one gain, or three global scalars with the simplified input pathway, leaving signs and routing discrete.Once the recurrence is structured, the ridge readout becomes the largest model object, so memory is governed by the task rather than reservoir size.Source: §4.2 · p. 14

5.1. Echo state property condition

• Exact orthogonality makes the recurrent contribution to the contraction bound exactly ρ_M, yielding a uniform-in-input ESP whenever C < 1.For the standard H-ESN, ρ_M < 1 is sufficient and ρ_M = 1 is the exact threshold under α = 1; the bound is sufficient but conservative for MF-H-ESN experiments.Source: §5.1 · pp. 16–17

5.2. Spectral and mixing properties

• At ρ = 0.99, the Hadamard operator achieves MC = 148 versus 151 for dense orthogonal and 45 for standard dense ESN, with the ordering holding over ρ ∈ [0.80, 0.99].The measurement uses tanh reservoirs at N = 256, with i.i.d. inputs, ridge readouts, and three seeds.Source: §5.2 · p. 19, Figure 5 · p. 21
• The Hadamard operator reaches PR = N = 256 after one step, whereas the SCR remains at PR = 1 because it shifts impulses without spreading them.The Hadamard and dense couplings later settle near PR ≈ N/3, while the SCR requires Θ(N) steps for full traversal.Source: §5.2 · p. 19, Figure 5 · p. 21
• The randomization elements prevent the bare Hadamard transform’s two degeneracies: involutive two-state cycling and cancellation of uniform drive into one coordinate.The permutation spreads eigenvalues around the circle, while sign diagonals prevent the uniform-drive collapse under coordinatewise nonlinearity.Source: §5.2 · p. 20

5. Theoretical analysis

• The Hadamard recurrence therefore matches dense orthogonal memory while providing dense-like one-step mixing with structured, multiplier-free operation.The comparison is made against SCR, dense orthogonal, and dense random recurrences at N = 256.Source: §5.2 · pp. 18, 20

5.3. Noise response

• The noise gain of orthogonal recurrences depends only on ρ and follows (1 − ρ^2)^−1/2 until saturation, while MF models are independent of ρ and dense gains vary by realization.The experiment injects σ = 10^−3 additive state noise at every step and measures stationary perturbed-clean distance normalized by single-step noise.Source: §5.3 · pp. 20–21

6.1. Experimental setup

• The experiments cover 20 classification and seven regression benchmarks, comparing Hadamard, dense, cycle, orthogonal, memristive-friendly, and GRU reservoirs under shared model-selection protocols.Reservoir hyper-parameters use 500 random-search configurations, ridge readouts, and three random initializations; datasets come from the UEA/UCR and Monash TSER archives.Source: §6.1 · pp. 22–23

6.2. Comparison at fixed reservoir size

• At N = 256, simplified-input Hadamard models achieve the two best mean classification accuracies, while Hadamard variants preserve dense-model performance across regression tasks.Simplified-input variants are best on 11 of 20 classification datasets; reservoir models outperform the GRU on the three FloodModeling regression tasks by factors of 2.5 to 13.Source: §6.2 · pp. 23–24

6.3. Scaling with the reservoir size

• From N = 256 to N = 8192, Hadamard models keep pace with dense reservoirs, while MF-H-ESNsi leads classification rankings and gains 2.2 accuracy points over SCR at N = 8192.Classification accuracy improves by 3.1–4.2 points for every model, whereas regression error remains nearly flat across sizes; MF-H-ESNsi is 1.7 points above dense MF-ESN at N = 8192.Source: §6.3 · pp. 27–28, Figure 7 · p. 28, Figure 8 · p. 29
• Three of four Hadamard variants rank among the top three reservoirs and are statistically indistinguishable from dense baselines across 27 benchmarks, confirming comparable accuracy at lower cost.The Friedman–Nemenyi analysis separates the GRU from the four best-ranked reservoirs, while reservoir models remain within one critical distance.Source: §6.3 · p. 29, Figure 9 · p. 30

6.4. Ablation of the structured operator

• The Hadamard transform supplies global mixing, permutation preserves memory near criticality, and sign diagonals preserve nonlinear state diversity under degenerate inputs.Without randomization, the bare transform has MC ≈2 and effective rank ≈2; at ρ_M = 0.95, removing permutation reduces memory from 198.1 ± 1.1 to 190.8 ± 4.4.Source: §6.4 · pp. 30–31, Figure 10 · p. 32

6.5. Robustness to quantization and noise

• On-device readout training recovers most state-noise losses, 8-bit quantization is nearly free, and permanent faults remain small: at most 2.3 points for 20% mismatch and 0.8 points for 5% sign flips.At σ = 1%, MF-H-ESN reaches 0.84 on short tasks versus 0.76–0.80 for leaky-tanh models; 5% routing faults cost less than half a point.Source: §6.5 · pp. 34–35, Table 3 · p. 35, Figure 12 · p. 36

Dense-coupling comparison

• The structured operator achieves whole-state one-step mixing with N log2 N additions, one programmed gain, and tens of kilobytes instead of dense N^2 analog values and 268 MB at N = 8192.Dense orthogonal recurrence preserves the same theoretical behavior but requires O(N3) factorization, dense storage, and dense application.Source: §7 · pp. 36–37

Overall comparison

• The best Hadamard variant ranks first among eight models on twenty classification benchmarks, while MF-ESN ranks first on regression and MF-H-ESN remains within a few percent of its error at every size.Table 4 summarizes the reference-model comparison at N = 8192; dense couplings trail the best Hadamard classification variant, while the memristive-friendly neuron drives regression performance.Source: Table 4 · p. 37, §7 · p. 37

Cycle comparison

• Compared with the cycle reservoir, the best Hadamard variants rank 3.4 and 3.7 versus 5.1 and 5.4 on classification and regression, respectively, while adding global mixing.The cycle uses N shifts per step and needs N steps to spread a coordinate impulse, whereas simplified-input Hadamard models use the same three global scalars as the cycle in analog values.Source: §7 · pp. 37–38

Implementation cost

• At N = 8192, the structured recurrence stores fifteen bits per unit and 15.4 kB overall, versus about 33 kB per unit for dense coupling.Its fixed Walsh-Hadamard transform stores nothing, but requires N log2 N additions per step and 2N log2 N physical connections.Source: §7 · p. 38

• O(N log N) operations, O(N) parameters, and no materialized matrix replace dense recurrent coupling while preserving representational properties and accuracy.The operator combines sign diagonals, a permutation, and a fast Walsh-Hadamard transform, removing the O(N^2) analog weights required by dense physical implementations.Source: §8 · p. 38
• Across twenty classification and seven regression benchmarks, Hadamard reservoirs track dense random orthogonal recurrences within one accuracy point at every size up to N = 8192 and increasingly outperform equally orthogonal cycle reservoirs.The critical-difference analysis separates only the fully trained GRU from the best-ranked reservoirs, supporting parity at a fraction of the cost.Source: §8 · pp. 38–39
• Global mixing explains the cycle-reservoir advantage, while the permutation protects memory, sign diagonals preserve nonlinear richness, and the bare involutive transform is degenerate.A coordinate impulse spreads across the whole state in one Hadamard step but never in the cycle; single-binary-wire variants rank best on classification and are most hardware-frugal.Source: §8 · p. 39
• 50× faster recurrent steps at N = 32768 and 15.4 kB versus 268 MB at N = 8192 quantify the structured operator’s device-dependent speed and device-independent memory advantages.The step overtakes dense multiplication from N ≈1024 on laptop CPUs and entry-level GPUs and from N ≈8192 on an A100; the simplified input pathway uses three analog values at every width.Source: §8 · p. 39
• Exact orthogonality makes the echo state condition tight, and simulations show essentially free eight-bit precision plus about two accuracy points of cost under 20% mismatch or discrete faults.The noise response is fixed at design time; flipped signs and rewired connections remain valid instances with the same spectrum and guarantees.Source: §8 · p. 39
• Open directions include physical memristive realization and aging, deep-layer composition, temporal parallel scans with structured projections, and identifying the linear Hadamard reservoir’s Walsh-domain representation.The within-step state mixing is composable with time-axis linear recurrences rather than an alternative, while the Walsh-domain representation remains unresolved.Source: §8 · p. 40

Implementation

• The dedicated CUDA kernel executes the full transform in one shared-memory kernel for N ≤8192, while larger sizes use a multi-stage global-memory implementation.It assigns one thread block per batch row and separates log2 N butterfly stages with block-level synchronization.Source: Appendix A · p. 41
• Pure TensorFlow launches log2 N kernels per recurrent step, whereas XLA fuses the butterfly stages into a small number of specialized GPU kernels.Both implementations use the radix-2 butterfly decomposition of the FWHT.Source: Appendix A · p. 41

Measurements

• Tables A.1–A.3 report platform-specific measurements summarized in Figure 3, with software and hardware versions documented in Appendix D.Source: Appendix A · p. 41

Accounting conventions

• The accounting separates stored parameters from trained parameters, counts recurrent, input, and bias pathways, and excludes the common readout from reservoir totals.Real-valued parameters use 32 bits, signs use one bit, and indices use log2 N or ⌈log2 d⌉ bits; global scalars are negligible in storage totals but counted as analog values.Source: Appendix B · pp. 41, 44

Operator implementation

• The FWHT requires N log2 N additions/subtractions and no multiplications, giving 4.9×10^4 additions versus 1.68×10^7 dense multiply-accumulates at N = 4096, a 341× reduction.Its butterfly network has log2 N stages with N/2 butterflies per stage and needs no coordinate reordering.Source: Figure A.1 · p. 42

Hardware timing

• On an Apple M3 Pro, operator-fused implementations cross over dense matmul at N ≈1024, while unfused FWHT crosses over at N ≈4096; on an A100, the crossover moves to N ≈8192.Below that A100 size, all implementations share a kernel-launch latency floor and XLA structured steps are never slower than dense ones.Source: Table A.1 · p. 43, Table A.3 · p. 43, Table A.2 · p. 43

Memory footprint

• At N = 8192, the simplified Hadamard reservoir occupies 17.4 kB versus 80.9 kB for H-ESN on a univariate task, while the operator accounts for 15.4 kB and 2N signs replace dense inputs.On MindReading (d = 204), the dense input pathway costs 6.7 MB versus 10 kB for the simplified pathway.Source: Appendix B · pp. 44–45, Table B.1 · p. 44

Readout size

• The five-class ridge readout uses (N + 1)C weights, occupying 20.5 kB at N = 1024 and 164 kB at N = 8192.At N = 8192, this readout is nine times the entire H-ESNsi reservoir and less than one-thousandth of the dense ESN.Source: Appendix B · p. 45

Trainable parameters

• At N = 256, GRUs optimize more parameters than reservoir readouts on all twenty benchmarks, with a median factor of 26 and up to 602 on FaceDetection.At N = 8192, the median ratio falls to 0.83, and GRUs are larger on nine of twenty benchmarks.Source: Appendix B · p. 45

• If γ ≤ 1 and εS+ ≤ γ, the state cube [0, H]N is forward-invariant under MF-H-ESN dynamics for any input.The componentwise update keeps each state between 0 and H.Source: Appendix C · p. 45
• Exact orthogonality makes the recurrent norm bound tight, yielding a uniform input-independent contraction with factor C; when C < 1, the ESP follows and trajectories converge at rate Ct.Orthogonality is the only operator-specific step that removes the slack in the contraction estimate.Source: Appendix C · p. 46
• The diagonal-free operator maps the constant input direction onto a single unit regardless of permutation, creating input degeneracy.This follows from the Hadamard transform’s constant first row and the signed diagonal construction.Source: Appendix C · pp. 46–47
• The full operator instead spreads constant-direction energy evenly across units in expectation, with norm-squared N, because the relevant sum contains N independent signs with zero mean and unit variance.Thus the full operator avoids the concentration produced by the diagonal-free variant.Source: Appendix C · p. 47

Benchmark scope

• The study uses 20 classification and 7 regression tasks spanning diverse classes, sequence lengths, input dimensionalities, dataset sizes, and sensor or biomedical recordings.Datasets were selected before inspecting results; impractical or too-small datasets were excluded.Source: Appendix D · p. 47

Implementation and timing

• Experiments use TensorFlow with aeon-loaded datasets, while model experiments run on an NVIDIA A100 and timing measurements use three platforms with warmed-up single-precision median executions.Timing advances 256 states per call; software versions are TensorFlow 2.12 on Apple M3 Pro and 2.21 on the NVIDIA T4 and A100.Source: Appendix D · p. 47

Reservoir search space

• Section 6.1 samples 500 configurations across recurrent scaling, input and bias gains, leak rates, memristive parameters, and ridge regularization, with physics constants fixed following.The searched ranges include ρ, ρM ∈ {0.7, 0.8, . . . , 1.2}, ωx, ωb, α, ε, γ, s, and ridge penalties from 10^-5 to 10^2.Source: Appendix D · p. 47

GRU training

• The GRU uses Adam, early stopping, and validation-based selection, while regression additionally expands the learning-rate grid because lower rates converged to the constant predictor.The wider regression-only grid includes 10^-2 and 3·10^-3; η = 10^-2 reduced training error to below one hundredth of the constant predictor’s error, whereas η ≤ 10^-3 did not escape the plateau.Source: Appendix D · pp. 47–48

Fixed-size comparison

• At fixed reservoir size, reservoirs use N = 256 while the GRU selects hidden size up to 256, and the protocol favors the trained baseline by retaining its validation split during assessment.Reservoir readouts are refit on the full training set after model selection, whereas GRU early stopping requires held-out data.Source: Appendix D · p. 48

Critical-difference analysis

• Critical-difference analysis selects reservoir size and GRU hidden size on the same validation split, allowing reservoirs up to N = 8192 while capping GRU hidden size at 256.Reservoir-size selection adds six values without further model selection, and the ranges remain comparable in trainable parameters.Source: Appendix D · p. 48

Classification benchmarks

• Classification curves show that accuracy saturates before N = 8192 on several benchmarks, degrades at the largest sizes on a few, and separates models more widely on harder tasks.Figures E.1 and E.2 show these twenty benchmark-specific curves; Figure 7 aggregates them.Source: Appendix E · p. 49

Regression benchmarks

• Regression curves show essentially flat error across N on most tasks, with Covid3Month spanning less than 6% of the error range, making apparent separations there meaningless.Figure E.3 reports the seven benchmark-specific curves.Source: Appendix E · p. 49