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

Andrea Ceni; Gianluca Milano; Carlo Ricciardi; Claudio Gallicchio

cs.AIAug 28, 2026

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

Authors: Andrea Ceni, Gianluca Milano, Carlo Ricciardi, Claudio Gallicchio

Organizations: Department of Computer Science, University of Pisa, Italy · Istituto Nazionale di Ricerca Metrologica (INRiM), Italy · Department of Applied Science and Technology, Politecnico di Torino, Italy

Abstract

Reservoir Computing (RC) designs Recurrent Neural Networks around a fixed, i.e., untrained, recurrent layer, and is a natural candidate for neuromorphic hardware. Memristive-friendly reservoirs derive the neuron dynamics from memristive-device kinetics, but still rely on dense recurrent matrices, which are expensive to realize physically. In this paper, we replace the dense matrix with a structured orthogonal operator, built from sign diagonals, a permutation, and a fast Walsh-Hadamard transform. The operator is multiplier-free, requires $O(N)$ parameters and $O(N\log N)$ operations per step, and is never materialized as a matrix. We instantiate it in a standard and in a memristive-friendly Echo State Network, with one binary input connection per unit. Our mathematical analysis shows that exact orthogonality yields an echo state condition that is tight in the recurrent scaling, and a noise response that is predictable at design time. Moreover, the operator mixes the whole state in a single application. Experiments on twenty classification and seven regression benchmarks, at reservoir sizes up to $N = 8192$ , show that the structured models match dense orthogonal reservoirs, and achieve better mean performance than the cycle reservoir by a margin that widens with size. Furthermore, we time the recurrent step on three hardware platforms, where it is up to $50\times$ faster than a dense product and $10^4\times$ smaller in memory. Finally, we ablate the operator and measure the response to noise, quantization, device mismatch and discrete faults.

Figures & tables

Figure 1: Comparison between the recurrent step of a conventional reservoir (a) and of the Hadamard reservoir proposed in this work (b). In (a) the dense matrix

\mathbf{W}_{h}

stores

N^{2}

analog weights and its application costs

N^{2}

multiply-accumulate operations. In (b) the same state mixing is obtained from the structured operator

\mathbf{M}=\mathbf{D}_{2}\widehat{\mathbf{H}}_{N}\Pi\mathbf{D}_{1}

(eq. 4 ), applied as a sequence of sign flips (

\mathbf{D}_{1}

\mathbf{D}_{2}

), a fixed routing (

\Pi

), and a fast Walsh-Hadamard transform (

\widehat{\mathbf{H}}_{N}

, detailed in A ). The stored description reduces to the

2N

signs of the two diagonals, the

N

indices of the permutation and the single scalar

\rho_{M}

, and the

N\log_{2}N

operations are additions and subtractions. No programmable multiplication is involved, as the recurrence is scaled by one constant gain, common to all the units, which can be absorbed into the entries of

\mathbf{D}_{2}

Figure 2: Computational graph of the memristive-friendly Hadamard reservoir (MF-H-ESN): the linear pre-activation of eq. 2 is replaced by the Hadamard version of eq. 6 , while the state update eq. 3 is unchanged. The dense recurrent term of the original memristive-friendly update [ 34 ] is replaced by the structured Hadamard operator

\rho_{M}\mathbf{M}

(compact box, detailed in Figure 1 b), which feeds the potentiation/depression kinetics

K_{p},K_{d}

and the leaky integration (dashed box) exactly as in the original model.

Figure 3: Wall-clock time per recurrent step (left) and speedup over the dense matmul baseline (right) as a function of the reservoir size

N

(fp32,

256

states advanced per call), on three hardware platforms. Top: laptop-class Apple M3 Pro CPU; middle: entry-level NVIDIA T4 GPU. On both, the structured recurrence crosses over the dense one at

N\!\approx\!1024

with a compiled implementation, and reaches about

50\times

N=32768

. Bottom: datacenter-class NVIDIA A100 GPU, where the crossover moves to

N\!\approx\!8192

. Below that size all implementations are dominated by a fixed per-call overhead, and the structured step stays within about ten percent of the dense one; beyond it the advantage grows (

3.5\times

N=16384

). Timings in tabular form are in A .

Figure 4: Memory footprint of the structured recurrence against the dense one, on a univariate task (

d=1

). Left: bits needed to store the recurrent operator and the input pathway, with real-valued parameters at

32

bits, signs at one bit and routing indices at

\log_{2}N

bits. Right: analog values that a physical realization has to program individually. The dense curve is the ESN, the dense orthogonal reservoir and the MF-ESN, which coincide, and each Hadamard curve is a

\tanh

and a memristive-friendly model, i.e., H-ESN and MF-H-ESN, and their simplified-input variants, which coincide in the same way, since they share the coupling and the input pathway; the analog values of the simplified-input variants are three at every size. The dotted line is the ridge readout of a five-class task, which every reservoir model carries in addition. Formulas, assumptions and the corresponding table are in B .

Figure 5: Memory and mixing at

N=256

(mean

\pm

std over 3 seeds). Left: memory capacity vs. the spectral radius

\rho

. The orthogonal recurrences (SCR, dense orthogonal, Hadamard) dominate the dense random ESN, and the Hadamard operator matches the dense orthogonal one. Right: participation ratio (eq. 10 ) of

\mathbf{M}^{t}\mathbf{e}_{i}

vs.

t

(log scale), at

\rho=0.95

and averaged over ten choices of the impulsed unit

i

;

\mathrm{PR}

does not depend on

\rho

. The SCR translates the impulse from one unit to the next without ever spreading it (

\mathrm{PR}=1

), whereas the Hadamard operator attains the maximum

\mathrm{PR}=N

in one step.

Figure 6: Empirical state-noise gain on synthetic signals, vs. the spectral radius

\rho

(

N=256

, mean

\pm

std over 3 seeds). The three orthogonal recurrences coincide, the dense random reservoir exceeds them at the edge of stability with a large seed variance, and the MF models show a flat gain, independent of

\rho

Figure 7: Scaling with the reservoir size

N=2^{8},\ldots,2^{13}

over the twenty classification benchmarks (3 instantiations each). Left: test accuracy averaged over datasets. Right: average rank among the eight models, computed per dataset (rank 1 is best; the axis is inverted so that better is up). Hyper-parameters are selected at

N=256

and reused across sizes for all models. The subscript si denotes the simplified-input variant defined in Section 3.2 . Per-dataset curves: E .

Figure 8: Scaling with the reservoir size on the seven regression benchmarks (3 instantiations each). Left: mean squared error averaged over benchmarks, logarithmic scale. Right: average rank among the eight models, computed per benchmark (rank 1 is best). Compare with Figure 7 : here the curves are nearly flat.

Figure 9: Critical difference diagram over the twenty-seven benchmarks (twenty classification, seven regression), with the reservoir size selected on validation for the reservoir models and the hidden size selected on validation for the GRU. Lower rank is better; models joined by a bar are not significantly different at the $5\%$ level (Friedman–Nemenyi, $\mathrm{CD}=2.31$ ; Friedman test on the nine models, $p=5\cdot 10^{-4}$ ).

Figure 10: The ablation variants as operators (

N=256

). Top: participation ratio (eq. 10 ) of

\mathbf{M}^{t}\mathbf{e}_{i}

\rho_{M}=0.95

, i.e., the number of active coordinates

t

steps after a coordinate impulse; the bare transform

\widehat{\mathbf{H}}_{N}

alternates between the whole state and a single unit, as its involutivity prescribes, and the permutation-only control

\Pi

never spreads the impulse. Bottom: memory capacity of the linear reservoir vs.

\rho_{M}

, with an inset zooming near the edge of stability, where removing the permutation lowers the memory and raises its variance, while removing the sign diagonals does not.

Figure 11: Sensitivity to a degenerate input map (

N=256

\rho_{M}=0.9

20

seeds). The single column of the input weight matrix

\mathbf{W}_{x}

is interpolated at fixed norm between i.i.d. entries (

\beta=0

) and the constant vector

\mathbf{1}

(

\beta=1

), along which every unit receives the same drive. Left: the effective rank of the resulting nonlinear states, i.e., the participation ratio (eq. 10 ) of their singular values, which counts the populated directions of the state space. It collapses for the operators without sign diagonals (

\widehat{\mathbf{H}}_{N}\Pi

and

\Pi

), while the full operator and

\mathbf{D}_{2}\widehat{\mathbf{H}}_{N}\mathbf{D}_{1}

(no permutation) keep it flat. Right: the memory capacity is flat for every variant that retains the Hadamard mixing, so the effect is purely nonlinear; only the mixing-free control

\Pi

degrades, and just at

\beta=1

. Cf. C .

	crossbar / butterfly mismatch $\sigma_{d}$				sign flips $p$			routing faults $p$
Model	$1\%$	$5\%$	$10\%$	$20\%$	$0.1\%$	$1\%$	$5\%$	$0.1\%$	$1\%$	$5\%$
ESN	$-0.1$	$+0.0$	$+0.0$	$-0.6$	—	—	—	—	—	—
SCR	$+0.1$	$+0.1$	$+0.5$	$+0.1$	—	—	—	—	—	—
Orth	$-0.0$	$+0.0$	$+0.2$	$-0.6$	—	—	—	—	—	—
MF-ESN	$-0.0$	$-0.0$	$-0.2$	$-0.4$	—	—	—	—	—	—
H-ESN	$-0.0$	$-0.2$	$-0.4$	$-2.3$	$-0.4$	$-0.4$	$-0.8$	$+0.0$	$-0.3$	$+0.3$

Table 3: Permanent device perturbations with on-device readout: change in test accuracy (percentage points, mean over the four benchmarks and

3

seeds) with respect to the unperturbed device. Crossbar mismatch applies

w\mapsto w(1+\delta)

to every entry of the materialized operator; butterfly mismatch applies it to the

2N\log_{2}N

nonzero entries of the FWHT stages; sign flips act on

\mathbf{D}_{1},\mathbf{D}_{2}

, and routing faults permute among themselves the destinations of a fraction of the units, leaving

\Pi

a permutation (both for the Hadamard models only).

Figure 12: Test accuracy under runtime perturbations at

N=256

(mean

\pm

standard error over 3 seeds; hyper-parameters from Section 6 , no re-tuning). Top row: mean over the three short-sequence tasks (ECG5000, Epilepsy, RacketSports), where the unperturbed accuracies lie between

0.91

and

0.93

; bottom row: FordA (

T=500

), where they lie between

0.61

and

0.71

and the dotted line marks the

0.5

level of this binary task. In the first two columns the perturbation grows from left to right, as the injected noise increases, in the deploy regime (a, d) and in the on-device regime (b, e); in the third column it shrinks from left to right, as the state is quantized to more bits (c, f). On the short tasks the on-device readout recovers most of the noise cost and the MF-H-ESN is the most robust model; on FordA the accumulated noise destroys the discriminative information for all models and recalibration no longer helps.

	ESN	SCR	Orth	MF-ESN	H-ESN,	H-ESN si ,
					MF-H-ESN	MF-H-ESN si
Operations per step	$N^{2}$ MAC	$N$ shift	$N^{2}$ MAC	$N^{2}$ MAC	$N\log_{2}N$ add	$N\log_{2}N$ add
Construction cost	$O(N^{2})$	$O(1)$	$O(N^{3})$	$O(N^{2})$	$O(1)$	$O(1)$
Reservoir memory, $N=8192$	$268$ MB	$2.1$ kB	$268$ MB	$268$ MB	$80.9$ kB	$17.4$ kB
Analog values in the reservoir	$N^{2}{+}2N$	$3$	$N^{2}{+}2N$	$N^{2}{+}2N$	$2N{+}1$	$3$
Impulse spread in one step	yes	no	yes	yes	yes	yes

Table 4: The proposed models against the reference ones, at a glance. The columns are, in order, the conventional dense ESN, the Simple Cycle Reservoir [ 35 ] , a dense random orthogonal recurrence, the memristive-friendly MF-ESN [ 34 ] , which shares the dense coupling of the ESN and differs in the neuron, and the models introduced here, in their two input variants; in the last two rows, where the two variants of a column differ, the entries refer to the

\tanh

and to the memristive-friendly neuron, in this order. The reservoir memory is the recurrent operator together with the input pathway and the bias, in the form in which they have to be resident to apply the recurrence, in single precision and on a univariate task; the analog values are the entries of the same three objects that a physical realization has to program one by one, and both are detailed in B . The readout, which is common to every reservoir model and occupies

164

kB in the same setting with five classes, is excluded from both. Mixing, echo state condition and noise gain are from Section 5 ; the last two rows are the mean rank among the eight models over the twenty classification and the seven regression benchmarks at

N=8192

, where rank

1

is the best (Section 6.3 ).

Appendix figures & tables8 assets

Supplementary material from the paper’s appendix.

Appendix

Figure A.1: Fast Walsh-Hadamard transform as a butterfly network, shown for

N=8

(

K=3

stages, stride

1,2,4

). Each butterfly combines two inputs into a sum (gray,

+

) and a difference (red,

-

). The coordinates need no reordering, neither before nor after the stages, so the network computes exactly the transform

\widehat{\mathbf{H}}_{N}

of Section 3.1 . In general, the

K=\log_{2}N

stages of

N/2

butterflies each amount to

N\log_{2}N

additions/subtractions and no multiplications, against

N^{2}

multiply-accumulates for the corresponding dense matrix-vector product: at

N=4096

4.9\times 10^{4}

additions against

1.68\times 10^{7}

multiply-accumulates, a

341\times

reduction.

$N$	dense [ms]	FWHT (TF)		FWHT (XLA)		FWHT (C++ op)
		[ms]	speedup	[ms]	speedup	[ms]	speedup
256	0.16	0.74	0.21 $\times$	0.31	0.51 $\times$	0.27	0.58 $\times$
512	0.38	1.16	0.32 $\times$	0.39	0.97 $\times$	0.42	0.90 $\times$
1024	1.19	2.23	0.53 $\times$	0.56	2.13 $\times$	0.72	1.64 $\times$
2048	3.28	4.30	0.76 $\times$	0.98	3.35 $\times$	1.47	2.23 $\times$
4096	11.43	10.09	1.13 $\times$	2.00	5.72 $\times$	2.85	4.01 $\times$

Appendix

Table A.1: Time per recurrent step [ms] and speedup vs. dense matmul (batch 256, fp32) on a laptop-class Apple M3 Pro CPU. The crossover occurs at

N\approx 1024

with operator fusion (XLA and C++ custom op) and at

N\approx 4096

for the unfused FWHT; the speedup profile closely mirrors that of the T4 GPU (Table A.2 ).

$N$	dense [ms]	FWHT (TF)		FWHT (XLA)		FWHT (CUDA op)
		[ms]	speedup	[ms]	speedup	[ms]	speedup
256	0.18	0.71	0.25 $\times$	0.34	0.53 $\times$	0.29	0.62 $\times$
512	0.38	1.18	0.32 $\times$	0.44	0.87 $\times$	0.43	0.88 $\times$
1024	1.16	2.12	0.55 $\times$	0.63	1.84 $\times$	0.74	1.57 $\times$
2048	3.28	4.16	0.79 $\times$	1.10	3.00 $\times$	1.57	2.09 $\times$
4096	11.74	10.50	1.12 $\times$	2.31	5.09 $\times$	3.08	3.81 $\times$

Appendix

Table A.2: Time per recurrent step [ms] and speedup vs. dense matmul (batch 256, fp32, NVIDIA T4).

$N$	dense [ms]	FWHT (TF)		FWHT (XLA)		FWHT (CUDA op)
		[ms]	speedup	[ms]	speedup	[ms]	speedup
256	0.20	0.52	0.39 $\times$	0.22	0.90 $\times$	0.30	0.66 $\times$
512	0.20	0.61	0.32 $\times$	0.23	0.87 $\times$	0.34	0.58 $\times$
1024	0.23	0.61	0.38 $\times$	0.24	0.98 $\times$	0.28	0.82 $\times$
2048	0.22	0.65	0.34 $\times$	0.23	0.94 $\times$	0.30	0.74 $\times$
4096	0.23	0.68	0.34 $\times$	0.24	0.97 $\times$	0.30	0.78 $\times$

Appendix

Table A.3: Time per recurrent step [ms] and speedup vs. dense matmul (batch 256, fp32, NVIDIA A100). On datacenter-class hardware the dense baseline is substantially faster than on the T4 (Table A.2 ), moving the crossover point from

N\approx 1024

N\approx 8192

; below that size all implementations sit on a common kernel-launch latency floor and the structured step is never slower than the dense one (XLA column).

	stored description (bits)		analog quantities		total, $d=1$
Model	recurrent	input, bias	values	connections	$N=1024$	$N=8192$
ESN	$32N^{2}$	$32N(d{+}1)$	$N^{2}{+}N(d{+}1)$	$N^{2}{+}Nd$	$4.2$ MB	$268.5$ MB
SCR	$32$	$N(d{+}1)$	$3$	$N{+}Nd$	$260$ B	$2.1$ kB
Orth	$32N^{2}$	$32N(d{+}1)$	$N^{2}{+}N(d{+}1)$	$N^{2}{+}Nd$	$4.2$ MB	$268.5$ MB
MF-ESN	$32N^{2}$	$32N(d{+}1)$	$N^{2}{+}N(d{+}1)$	$N^{2}{+}Nd$	$4.2$ MB	$268.5$ MB
H-ESN	$N(2{+}\log_{2}N)$	$32N(d{+}1)$	$1{+}N(d{+}1)$	$2N\log_{2}N{+}Nd$	$9.7$ kB	$80.9$ kB

Appendix

Table B.1: Memory footprint of the reservoir description. Left: bits required to store the recurrent operator and the input pathway, with real-valued parameters at

32

bits, signs at one bit, and indices at

\log_{2}N

and

\lceil\log_{2}d\rceil

bits; the few global scalars per model are negligible here and are omitted, but they are counted among the analog values, which they exhaust for the simplified-input models. Middle: quantities that a physical realization has to instantiate, i.e., analog values to be programmed individually, and analog connections to carry them; the butterfly interconnect of the FWHT and the cycle of the SCR are wiring, and are counted here only. Right: resulting total for the reservoir description on a univariate task (

d=1

), at two reservoir sizes. The readout, which is common to all the reservoir models, is listed separately for a

C

-class task, with

C=5

Figure E.1: Test accuracy (mean

\pm

std over 3 instantiations) versus reservoir size, one panel per classification benchmark (first ten, in alphabetical order; the remaining ten in Figure E.2 ).

Figure E.2: Test accuracy versus reservoir size, remaining ten classification benchmarks. Continuation of Figure E.1 .

Figure E.3: Test mean squared error (mean

\pm

std over 3 instantiations, logarithmic scale) versus reservoir size, one panel per regression benchmark. Lower is better. Note the different vertical ranges: on Covid3Month all models lie within a few percent of each other, as discussed in Section 6.2 .

Reservoir Computing

Reservoir computing uses recurrent dynamical systems to transform temporal inputs into features that support efficient learning with simple readouts.

See all 107 papers

On the Role of Preprocessing and Memristor Dynamics in Reservoir Computing for Image Classification
Apr 23, 2026Rishona Daniels, Duna Wattad, Ronny Ronen +2Reservoir Computing Neuromorphic Computing
Viterbi Faculty of Electrical and Computing Engineering, Technion - Israel Institute of Technology, Haifa, Israel · College of Engineering and Physical Sciences, Aston University, Birmingham B4 7ET, United Kingdom
Lindblad-Inspired Multi-Timescale Reservoir Computing with Separable Rotation and Dissipation
Jul 27, 2026Jyotiranjan Beuria, Amit ShuklaDynamic Neural Networks Reservoir Computing
IKS Research Centre, ISS, Delhi, India · Indian Institute of Technology Mandi, Himachal Pradesh, India
Scalable Memristive-Friendly Reservoir Computing for Time Series Classification
Apr 21, 2026Coşku Can Horuz, Andrea Ceni, Claudio Gallicchio +1Dynamic Neural Networks Energy-Efficient ML
Adaptive AI Research Group University of Lübeck, Germany · Department of Computer Science University of Pisa, Italy

Apr 23, 2026J/K move · Enter open · S save

On the Role of Preprocessing and Memristor Dynamics in Reservoir Computing for Image Classification

Rishona Daniels, Duna Wattad, Ronny Ronen, David Saad+1

Viterbi Faculty of Electrical and Computing Engineering, Technion - Israel Institute of Technology, Haifa, Israel · College of Engineering and Physical Sciences, Aston University, Birmingham B4 7ET, United Kingdom

Reservoir Computing Neuromorphic Computing Image Classification

Open paper AskPDF arXiv

Reservoir computing (RC) is an emerging recurrent neural network architecture that has attracted growing attention for its low training cost and modest hardware requirements. Memristor-based circuits are particularly promising for RC, as their intrinsic dynamics can reduce network size and parameter overhead in tasks such as time-series prediction and image recognition. Although RC has been demonstrated with several memristive devices, a comprehensive evaluation of device-level requirements remains limited. In this paper, we analyze and explain the operation of a parallel delayed feedback network (PDFN) RC architecture with volatile memristors, focusing on how device characteristics -- such as decay rate, quantization, and variability -- affect reservoir performance. We further discuss strategies to improve data representation in the reservoir using preprocessing methods and suggest potential improvements. The proposed approach achieves 95.89% classification accuracy on MNIST, comparable with the best reported memristor-based RC implementations. Furthermore, the method maintains high robustness under 20% device variability, achieving an accuracy of up to 94.2%. These results demonstrate that volatile memristors can support reliable spatio-temporal information processing and reinforce their potential as key building blocks for compact, high-speed, and energy-efficient neuromorphic computing systems.

Abstract

Reservoir Computing

Explore similar work

On the Role of Preprocessing and Memristor Dynamics in Reservoir Computing for Image Classification

Lindblad-Inspired Multi-Timescale Reservoir Computing with Separable Rotation and Dissipation

Scalable Memristive-Friendly Reservoir Computing for Time Series Classification