Free-Probability Kernels for Zero-Rollout Hyperparameter Selection in Reservoir Computing

1

Free-Probability Kernels for Zero-Rollout Hyperparameter Selection in Reservoir Computing

arXiv:2608.20998v1 [cs.LG] 21 Aug 2026

Sara Malacarne1 , Andrea Ceni2 , Claudio Gallicchio2 1 Telenor Research & Innovation, Oslo, Norway 2 Department of Computer Science, University of Pisa, Pisa, Italy

Abstract—Reservoir computing (RC) couples a fixed recurrent dynamical system with a trained lightweight readout, but this efficiency is partly lost during hyperparameter selection: the recurrent gain, input scale, and leakage rate determine the reservoir’s stability and temporal processing regime and are usually tuned through many rollouts. We introduce a deterministic, pilot-informed selector for leaky linear reservoirs followed by coordinate-wise nonlinear features. Free probability yields crosslag propagation coefficients that summarize how the reservoir mixes past inputs. In the large-width limit, these coefficients define a deterministic temporal kernel that approximates the finite-reservoir feature geometry. Kernel ridge regression on a short labeled pilot sequence therefore ranks candidate operating regimes without instantiating or rolling out a reservoir, and the selected configuration transfers across widths. Across ten synthetic temporal benchmarks, zero-rollout selection obtains a mean deployment score of 0.772, compared with 0.774 for exhaustive simulation-based search, while avoiding 156 600 selection rollouts. With a small rollout budget, the proposed ranking provides the strongest mean performance at every tested budget and reaches the exhaustive reference using 4.8% of its rollout cost. On four public electricity-transformer-temperature (ETT) forecasting datasets, five retained candidates recover the exhaustive operating point on three datasets. On multivariate cellular-traffic forecasting, 15 rollouts per cell reach the 462rollout exhaustive reference and outperform random search and Bayesian optimization at low budgets. These results position free-probability kernels as deterministic surrogates for selecting reservoir operating regimes when validation rollouts are scarce. Index Terms—Reservoir computing, Echo state networks, Free probability, Kernel methods, Hyperparameter selection

I. I NTRODUCTION Reservoir computing (RC) provides an efficient framework for temporal sequence learning: a fixed recurrent system maps the input history into a high-dimensional representation, while only a lightweight readout is trained [1], [2]. Its canonical instance is the Echo State Network (ESN) [3]. Although readout training is inexpensive, performance depends strongly on hyperparameters controlling the recurrent dynamics, particularly the recurrent gain, input scale, and leakage rate. These parameters govern stability, memory, and the relative influence of recent and distant inputs, and their best values can vary substantially across tasks and temporal scales. In standard practice, these hyperparameters are selected by repeatedly instantiating finite-width reservoirs, generating their state trajectories, fitting readouts, and evaluating validation error. Consequently, the search cost grows with the candidate grid, the reservoir width, and the number of random realizations. The selected configuration may also depend on the width

used during search, so selection may need to be repeated when the deployment width changes. We address this bottleneck with a deterministic, taskinformed selector that requires no finite-reservoir rollout during candidate ranking. Given a short labelled input–output sequence from the downstream task, the method evaluates each candidate through a deterministic large-width kernel. The pilot sequence describes the prediction task; it is not obtained by running the candidate reservoirs. The resulting ranking can be used directly as a zero-rollout selector or as a pre-screen that reduces the candidate grid before a small number of task-direct finite-reservoir evaluations. We study a leaky linear recurrence followed by coordinatewise nonlinear readout features. For this architecture, the nonlinear feature Gram matrix converges at large reservoir width to a deterministic temporal kernel. Its covariance is governed by mixed cross-lag propagation moments of the recurrent matrix, which we compute using free probability (FP). For each operating point θ = (σr , σin , α), the resulting kernel approximates the feature geometry of the finite reservoir with those same parameters. Kernel ridge regression on the pilot sequence can therefore rank the candidate operating points, after which the selected parameters are transferred to the finite reservoir used for deployment. Existing reservoir kernels are used mainly either to describe large-width reservoir behaviour or as predictors in their own right. Our use is different: for each candidate θ = (σr , σin , α), the deterministic kernel approximates the feature Gram matrix of the finite reservoir with the same parameters. We use this approximation to rank candidates on labelled pilot data, and then deploy the selected θ⋆ in the finite reservoir. The main analytical difficulty is that the same recurrent matrix is applied repeatedly over time, so the effects of inputs arriving at different lags are not independent. Our derivation captures these cross-time dependencies and shows that, at large width, the resulting nonlinear feature similarities converge to a deterministic kernel. To the best of our knowledge, this is the first task-informed, zero-rollout method that uses such a kernel to select the hyperparameters of a finite reservoir. The main contributions are as follows: 1) We derive the deterministic large-width post-nonlinearity kernel of a leaky linear reservoir. The analysis computes the mixed cross-lag propagation moments generated by repeated use of the same recurrent matrix and establishes the coordinate-level self-averaging required by the nonlinear feature map. We also obtain a controlled complete-

2

history extension and explicit kernels for several structured recurrent ensembles. 2) We use this kernel as a task-informed surrogate for finitereservoir hyperparameter selection. Candidate values of (σr , σin , α) are ranked from labelled pilot data without instantiating any candidate reservoir. If the deterministic validation score has a unique best candidate, then a sufficiently wide finite reservoir selects the same candidate with high probability. 3) We evaluate zero-rollout selection and FP-based prescreening on synthetic temporal benchmarks, the public ETT forecasting datasets [4], and an operational cellulartraffic forecasting task. We compare against exhaustive task-direct selection, a memory-based proxy, nonlinearESN selection, random search [5], and Bayesian optimization using the tree-structured Parzen estimator (TPE) [6], under matched finite-reservoir rollout budgets. Complete proofs, an explicit kernel catalogue for alternative recurrent random matrices, and additional experimental results are provided in Appendices A–G. II. BACKGROUND AND P ROBLEM F ORMULATION Existing approaches to reservoir hyperparameter selection can be divided broadly into task-based search and dynamical proxies. Task-based methods evaluate candidate configurations according to their downstream validation performance. The simplest examples are grid and random search [5], while more adaptive strategies include Bayesian optimization for reservoir parameters [7], gradient-based optimization [8], and evolutionary methods [9]. More efficient validation schemes reduce or reuse some of these evaluations [10], [11]. Despite their different search strategies, these methods share the same basic cost: to assess a new candidate, they instantiate a finite reservoir, generate its state trajectory, fit a readout, and evaluate its predictions. A different strategy is to select or adapt reservoirs using properties of their internal dynamics. Echo-state-property conditions and spectral-radius guidelines identify stable operating regions [3], [12], [13], while memory capacity and Jacobianbased indicators quantify aspects of temporal processing [14]. Reservoir dynamics may also be adapted without labels, for example through intrinsic plasticity [15]. These criteria provide useful information about stability, memory, or dynamical regime, but they are not designed to rank candidates directly according to the labelled downstream validation objective. When they depend on realized states, they also still require a finite reservoir to be instantiated and run. Large-width reservoir theory offers a third perspective: replacing the random finite system by a deterministic limiting description. Hermans and Schrauwen formulated infinite-width ESNs as recurrent kernel machines [16]. Couillet et al. used random-matrix theory to characterize the training and testing behaviour of large linear ESNs [17], while Dong et al. derived nonlinear recurrent kernels and used them directly for prediction [18]. Related recurrent-kernel constructions have since been developed for leaky, sparse, and deep reservoir topologies [19]. Gonon, Grigoryeva, and Ortega [20] instead

construct an infinite-dimensional Volterra reservoir whose kernel defines a universal temporal predictive model. These works show that the behaviour of a large reservoir can often be represented by a deterministic kernel. Their main purpose, however, is to characterize the reservoir or to use the limiting kernel itself as the predictive model. We address a different question: whether such a kernel can replace finite-reservoir rollouts during hyperparameter selection. For each candidate θ = (σr , σin , α), we construct the large-width kernel corresponding to the finite reservoir with those same parameters. Its performance on labelled pilot data is used to rank the candidates, after which the selected θ⋆ is transferred back to the finite reservoir used for deployment. We now introduce the reservoir family for which this correspondence is derived. Consider the leaky linear recurrence xt+1 = Axt + αWin ut ,

A = (1 − α)In + αWr ,

(1)

equipped with coordinate-wise nonlinear readout features zt = ψ(xt ),

⊤ ŷt = Wout zt .

(2)

Here, ut ∈ Rdin , xt , zt ∈ Rn , and ŷt , yt ∈ Rdout denote the input, linear state, nonlinear feature vector, prediction, and target, respectively. The recurrent and input matrices Wr ∈ Rn×n and Win ∈ Rn×din are randomly initialized and remain fixed, while only Wout ∈ Rn×dout is fitted. The leakage rate is α ∈ (0, 1], the reservoir is initialized at x0 = 0, and ψ is applied coordinate-wise. The experiments use ψ = tanh. The scalings of Wr and Win introduce the recurrent gain σr and input scale σin and are specified in Section III. Placing the nonlinearity after the recurrence separates temporal memory from nonlinear processing. Related lineartransition designs also appear in modern state-space sequence models such as S4 and LRU [21], [22]. Linear recurrent dynamics do not by themselves imply limited temporal modelling capacity: stable linear reservoirs equipped with sufficiently expressive polynomial or neural-network readouts form universal approximation families for broad classes of fading-memory filters [23], [24]. In such architectures, the linear recurrence represents and propagates the input history, while the readout supplies the nonlinear approximation capacity. Our model uses the lightweight feature map zt = tanh(xt ) followed by a trained linear output layer. This coordinatewise readout is not itself covered verbatim by the broader universality theorem, but it has recently been shown to be competitive with fully nonlinear ESNs on a range of temporal tasks [25]. More importantly for the present analysis, placing the nonlinearity outside the recurrent loop keeps the state linear in the propagated input history. Indeed, unrolling (1) from x0 = 0 gives xt = α

t−1 X

Ak Win ut−1−k .

(3)

k=0

Because the same recurrent matrix is applied repeatedly over time, the covariance between states depends on the mixed normalized moments
1 Tr Ak (Aℓ )⊤ . (4) n

3

These moments form the reservoir’s cross-lag propagation profile: diagonal terms measure the average squared gain of delayed echoes, whereas off-diagonal terms measure their overlap. Unlike the spectral radius, this profile also retains cross-time coupling and non-normal amplification. The next section derives the large-width limits of these mixed moments, uses them to construct the nonlinear feature kernel, and applies the kernel to zero-rollout candidate ranking. III. F REE -P ROBABILITY K ERNEL FOR T EMPORAL M ODEL S ELECTION Sections III-A and III-B specify the reservoir family and introduce the empirical-kernel notation. The new theoretical results begin in Section III-C, where we derive the deterministic large-width kernel at fixed context length; Sections III-D and III-E extend this construction to the complete washed-out history and use it for hyperparameter selection. Throughout this section, θ = (σr , σin , α) denotes a fixed candidate configuration. To avoid overloading the notation, we suppress the dependence on θ in A, Wr , Win , and the reservoir states, but retain it in the kernel and covariance quantities because the selector compares these objects across candidates. The subscript n identifies a finite-width random kernel, while kernels without n are deterministic large-width limits. The superscript (L) denotes an L-context quantity, and (∞) denotes complete history.

by the same input sequence but initialized at x0 and x′0 , respectively. Their difference satisfies xt − x′t = At (x0 − x′0 ). Hence, the influence of the initial condition vanishes for every pair x0 , x′0 if and only if At → 0 as t → ∞, which is equivalent to ρ(A) < 1. (8) √ For Wr = (σr / n)Gr , the circular law places the largewidth eigenvalue cloud of Wr in a disk centered at the origin with radius σr [26]. Therefore, the eigenvalue cloud of A = (1 − α)In + αWr is asymptotically a disk centered at 1 − α with radius ασr . Its asymptotic spectral radius is consequently ρFP (α, σr ) := (1 − α) + ασr .

For α > 0, the asymptotic linear echo-state condition ρFP < 1 is equivalent to σr < 1. A finite-width robustness margin is introduced in Section III-E. B. Readout and finite-reservoir kernel notation The matrices Wr and Win in eq. (5) are never trained. The readout features are zt = ψ(xt ),

A. Leaky linear reservoir and stability We use the leaky linear reservoir model of Section II, eq. (1) with a scaled real Ginibre recurrent matrix Wr and an independent Gaussian input matrix Win : σr Gr ∈ Rn×n , Wr = √ Gr , n σin Gin , Gin ∈ Rn×din , Win = √ (5) din i.i.d.

(Gr )ij ∼ N (0, 1),

i.i.d.

(Gin )ij ∼ N (0, 1),

Gr ⊥ Gin . √

The 1/ n normalization keeps the spectrum and operator norm of Wr on an O(1) scale as the reservoir width grows. The normalization of Win controls the scale of the injected input. Indeed, for fixed deterministic vectors u, v ∈ Rdin , 2   σin ⊤ EWin Win uv ⊤ Win = (u⊤ v)In . din

(6)

In particular, VarWin [(Win u)i ] =

2 σin ∥u∥2 . din

(7)

Our asymptotic analysis takes n → ∞ with din fixed; no additional high-dimensional assumption on the input vectors is required. The independence of Gr and Gin allows us to condition on the recurrent matrix while retaining an independent Gaussian input map. For the linear recurrence, the echo-state property can be characterized exactly. Let xt and x′t be two trajectories driven

(9)

zt,i = ψ(xt,i ),

i = 1, . . . , n,

(10)

where ψ : R → R is applied coordinate-wise. The experiments use ψ = tanh. The theoretical kernel result assumes ψ ∈ Cb1 (R), meaning that ψ is bounded and continuously differentiable with bounded derivative. Let T and V denote the training and validation time indices of the pilot sequence. Stacking the corresponding untruncated feature vectors gives

ZT = zt⊤ t∈T , ZV = zt⊤ t∈V . (11) Since the feature map ψ is fixed within each construction, we suppress it in the notation for the pilot kernel blocks. The width-normalized empirical kernel blocks for candidate θ are 1 ZT ZT⊤ , n 1 Kn,θ,V T = ZV ZT⊤ . n Kn,θ,T T =

(12)

They are indexed by pilot time steps, not by reservoir coordinates. Writing YT for the training targets, kernel ridge regression predicts on the validation block as −1 YbV,n (θ, λ) = Kn,θ,V T (Kn,θ,T T + λI) YT .

(13)

The blocks in (12) are computed from the untruncated finite-reservoir trajectory. Section III-C first derives deterministic large-width limits for their L-context counterparts. Section III-D then shows that, in the stable regime, the finite-context kernels converge to the complete-history kernel associated with the washed-out recurrence.

4

C. Finite-context free-probabilistic kernel construction Starting from x0 = 0, the linear recurrence can be unrolled exactly: t−1 X xt = α Ak Win ut−1−k . (14) k=0

For a fixed context length L < ∞, the selector retains the most recent L + 1 propagated inputs and uses the truncated state L X (L) xt := α Ak Win ut−1−k . (15) k=0

The corresponding empirical nonlinear kernel entry at width n is 1 (L) (L) (16) Kψ,n,θ (t, s) := ψ(xt )⊤ ψ(x(L) s ). n For the pilot train and validation sets, define the corresponding empirical L-context blocks entrywise by  (L)  (L) Kn,θ,T T t,s := Kψ,n,θ (t, s), t, s ∈ T , (17)  (L)  (L) Kn,θ,V T t,s := Kψ,n,θ (t, s), t ∈ V, s ∈ T . Because the recurrence is linear, the state covariance is determined by mixed normalized traces of propagated recurrent matrices. For fixed k, ℓ ≥ 0, define min(k,ℓ) 

τk,ℓ :=

X j=0

k j

  ℓ (1 − α)k+ℓ−2j (ασr )2j . j

(18)

These coefficients quantify how input components injected at two lags remain correlated after propagation through the reservoir. Proposition 1 (Ginibre cross-lag propagation coefficients). Let σr An = (1 − α)In + αWr(n) , Wr(n) = √ G(n) . n r For every fixed k, ℓ ≥ 0,
L1 1 Tr Akn (Aℓn )⊤ −−→ τk,ℓ . (19) n The proof expands the mixed powers by the binomial theorem, uses the circular-element limit for the normalized ∗moments, and controls fluctuations by Gaussian concentration. Complete details are given in Appendix A. The trace limit controls only an average over coordinates. Because the activation is applied coordinate-wise, we additionally require the diagonal propagation terms to self-average. Lemma 1 (Averaged diagonal self-averaging). Let An = (1 − (n) α)In + αWr . For every fixed L < ∞, n

2 P 1 X  k ℓ ⊤ An (An ) ii − τk,ℓ − → 0. max 0≤k,ℓ≤L n i=1

(20)

Thus, although individual reservoir units remain random, their diagonal propagation profile becomes deterministic. This coordinate-level result is essential because ψ acts neuronwise; dependence among different units is handled by the concentration step below. The complete proof is given in Appendix A.

Using the cross-lag propagation coefficients, define the deterministic linear-state covariance L

(L)

Qθ (t, s) := α2

2 X σin τk,ℓ u⊤ t−1−k us−1−ℓ , din

(21)

k,ℓ=0

where θ = (σr , σin , α), and the associated 2 × 2 covariance matrix ! (L) (L) Qθ (t, t) Qθ (t, s) (L) Σθ (t, s) := . (22) (L) (L) Qθ (t, s) Qθ (s, s) The nonlinear feature kernel is obtained by passing the limiting Gaussian state covariance through the coordinate-wise activation. Theorem 1 (Deterministic nonlinear reservoir kernel). Let ψ ∈ Cb1 (R). For fixed L, t, s and a deterministic bounded input sequence (ur )r , as n → ∞, (L)

P

(L)

Kψ,n,θ (t, s) − → Kψ,θ (t, s), where

(23)

(L)

Kψ,θ (t, s) := E[ψ(gt )ψ(gs )] ,   (L) (gt , gs ) ∼ N 0, Σθ (t, s) .

(24)

The convergence is in probability with respect to the joint randomness of Wr and Win . a) RC interpretation and proof idea: Conditional on the recurrent matrix, each reservoir unit receives a Gaussian projection of the same propagated input history. Lemma 1 shows that the unitwise covariance profiles self-average, and concentration over the input weights removes the remaining finite-width fluctuations. Hence a wide random reservoir induces a deterministic temporal feature geometry at kernel level. On any fixed pilot set, the empirical Gram blocks converge in Frobenius norm to their deterministic counterparts. Complete details are given in Appendix A. The nonlinear expectation in (24) can be evaluated for the exact feature map ψ = tanh by deterministic lowdimensional Gaussian quadrature. For each pair (t, s), the (L) covariance matrix Σθ (t, s) defines a bivariate Gaussian law, (L) and Ktanh,θ (t, s) is obtained by Gauss–Hermite quadrature of E[tanh(gt ) tanh(gs )]. Since Theorem 1 holds for any ψ ∈ Cb1 , a closed-form erf/arcsine surrogate is also available; its uniform entrywise approximation bound and empirical validation are given in Remark 1 and Appendix F-A. The experimental use of the exact kernel and of this surrogate is specified in Section IV. D. Complete history and other recurrent ensembles A finite context is the implementable selector; the completehistory limit connects it to the fading-memory interpretation of RC. When σr < 1, the reservoir forgets remote inputs geometrically, so increasing L recovers the washed-out state with controlled error. For a bounded input history and σr < 1, define ∞ σ2 X (∞) τk,ℓ u⊤ (25) Qθ (t, s) := α2 in t−1−k us−1−ℓ , din k,ℓ=0

5

(∞)

and let Kψ,θ (t, s) be the Gaussian nonlinear kernel obtained from this covariance as in (24). The covariance series converges absolutely, and the finite-context kernels converge geometrically to the complete-history kernel: (L)

(∞)

Kψ,θ (t, s) −→ Kψ,θ (t, s),

L → ∞.

(26)

Because the admissible candidate grid is finite, this convergence is uniform over the candidates used by the selector. For a complete bounded input history, let (∞)

xt

:= α

∞ X

Ak Win ut−1−k

(27)

k=0

denote the complete washed-out state. For real Ginibre recurrence with σr < 1, this series is well defined almost surely for all sufficiently large n. Strong convergence controls the powers of the non-normal transition matrix and hence the finite-width truncation tail. Consequently, on every fixed finite pilot set, 1 P (∞) (∞) ) −−−−→ Kψ,θ (t, s). ψ(xt )⊤ ψ(x(∞) s n→∞ n (28) The quantitative truncation bounds and the full Ginibre argument are given in Appendix B. The same construction applies to recurrent ensembles with deterministic mixed propagation moments and diagonal selfaveraging. The Haar-orthogonal argument and explicit kernels for cyclic, circulant, skew-symmetric, and complex-valued recurrences are given in Appendices C–D. The experiments in this paper use real Ginibre recurrence throughout. (∞)

Kψ,n,θ (t, s) :=

E. Zero-rollout selection on temporal pilot data The kernel now acts as a virtual reservoir during model selection. For each operating point, it captures the large-width temporal feature geometry of the corresponding finite reservoir without constructing recurrent weights or state trajectories. For each candidate θ = (σr , σin , α), the construction above replaces the empirical feature kernel of a large random reservoir by a deterministic kernel computed directly from a short (L) labelled pilot sequence. Let Ksel,θ denote the deterministic selection kernel used in a given experiment. In the experiments below, this is either the exact tanh kernel evaluated by quadrature or the closed-form erf surrogate described above; the precise choice is stated in Section IV. For pilot train and validation index sets T and V, define the deterministic kernel blocks entrywise by  (L)  (L) Kθ,T T t,s := Ksel,θ (t, s), t, s ∈ T , (29)  (L)  (L) Kθ,V T t,s := Ksel,θ (t, s), t ∈ V, s ∈ T , (L)

where each entry is computed from Qθ by either exact tanh quadrature or the closed-form erf surrogate, as specified above. For a ridge parameter λ, the corresponding validation prediction is  −1 (L) (L) YbV (θ, λ) = Kθ,V T Kθ,T T + λI YT . (30)

Algorithm 1: Zero-rollout FP hyperparameter selector Input: Pilot sequence (ut , yt ), split into train/validation; context length L; candidate grid Θ; ridge grid Λ Output: Selected hyperparameters θ⋆ and ranking regularizer λ⋆sel for each θ = (σr , σin , α) ∈ Θ satisfying the stability guard do (L) 2 Compute τk,ℓ (θ) via (18) for 0 ≤ k, ℓ ≤ L; Build Qθ (L) (L) (L) via (21); compute Ksel,θ ; form Kθ,T T and Kθ,V T via (29); for each λ ∈ Λ do 3 Compute YbV (θ, λ) via (30); compute (L) NMSEval (θ, λ); 4 end 5 end (L) ⋆ ⋆ 6 return (θ , λsel ) = arg minθ,λ NMSEval (θ, λ); 1

a) FP selector: Algorithm 1 presents the holdout-ridge version used for the synthetic benchmarks. For ETT and Telco, the inner selection of λ is instead performed by generalized cross-validation on the training block, after which the candidate is scored on the chronological validation block. The reservoir-candidate ranking procedure is otherwise unchanged. All matrices in Algorithm 1 are indexed by pilot time steps rather than reservoir coordinates: no n × n recurrent matrix is formed and P no finite reservoir P is rolled out. Here (L) NMSEval (θ, λ) = t∈V ∥b yt − yt ∥2 / t∈V ∥yt − ȳV ∥2 is the normalized mean-squared error on the pilot validation block, with ȳV the validation target mean. After precomputing input 2 L2 ) inner products, each candidate kernel requires O(Tpilot 2 lag-pair operations and O(Tpilot ) memory, followed by kernel ridge regression on the Ttrain × Ttrain block. Crucially, none of these costs depends on the deployment width n. The regularizer λ⋆sel is used only for ranking; at deployment, λ is re-selected on the deployment training block over the same canonical grid, using the kernel-scale feature convention of Section IV, so that the selection-stage and deployment-stage regularizers are directly comparable. Theorem 1 concerns kernel entries; the selector, however, returns an arg min over a finite grid. The following corollary records the corresponding decision-level result. Corollary 1 (Consistency of FP selection). Assume that the deterministic kernel and the finite reservoir use the same feature map. Fix finite candidate and ridge grids with strictly positive ridge parameters. For any fixed context length L, if the deterministic validation score has a unique minimizer, then finite-width reservoir selection returns the same minimizer with probability tending to one as n → ∞. In the stable regime, a unique complete-history optimum is preserved by all sufficiently large finite contexts, and selection based on the complete finite-reservoir trajectory converges to this optimum as n → ∞. The result follows from kernel convergence, continuity of kernel ridge regression, and stability of a unique minimizer on a finite grid. The complete statement and proof are given in Appendices A–B. The corollary concerns agreement on the pilot validation

6

objective; it does not assert that the selected configuration minimizes deployment test error. When several candidates are nearly tied, agreement may require a larger context length or reservoir width. In the synthetic experiments, the selector uses the erf surrogate while the finite reservoir uses tanh, so the corollary does not apply verbatim. The selection-level impact of this approximation is evaluated separately in Remark 1 and Appendix F-A. b) Stability guard: Equation (9) gives the asymptotic spectral radius of the Ginibre transition matrix. Since the largewidth selector cannot detect realization-specific instabilities of finite Ginibre matrices close to the stability boundary, we apply the fixed conservative guard ρFP (α, σr ) ≤ 1 − δ,

δ = 0.05.

(31)

The margin is a practical finite-width safeguard rather than a claim of optimality, and it is kept fixed across deployment widths. We apply the guard to every candidate grid before selection. For recurrent matrices whose finite-width spectral boundary is controlled exactly, such as scaled Haar-orthogonal or cyclic-shift recurrence, the additional margin is not required. IV. E XPERIMENTS A. Setup We evaluate the selectors in three complementary regimes. First, a suite of ten synthetic sequence benchmarks covers nonlinear temporal processing, explicit memory, delayed nonlinear transformations, and forecasting of chaotic dynamical systems. Second, four public ETT-small datasets (ETTh1, ETTh2, ETTm1, and ETTm2) test open-data real-world forecasting at hourly and 15-minute sampling rates. Third, a proprietary cellular-traffic dataset tests operational multivariate forecasting on ten randomly selected network cells. Thus, the evaluation spans controlled mechanistic tasks, reproducible public forecasting data, and a cell-replicated operational application. Synthetic benchmarks: The synthetic suite contains NARMA-k for k ∈ {10, 20, 30, 50} [27], linear memory capacity (MC) [28], the nonlinear delayed-input task of Inubushi and Yoshimura [29], and four forecasting tasks generated from chaotic dynamical systems: Lorenz63 [30], Mackey–Glass with two delay/horizon settings [31], and 5D Lorenz96 [32]. Their deployment protocols and task parameters are summarised in Table I. a) Real-data benchmarks: We use the public ETT-small benchmark introduced in [4], comprising ETTh1, ETTh2, ETTm1, and ETTm2. For ETT, the input consists of oil temperature (OT) and calendar covariates, and the target is future OT at physical horizons of 1, 6, and 12 hours. ETTh1/ETTh2 are sampled hourly, whereas ETTm1/ETTm2 are sampled every 15 minutes. The proprietary Telco benchmark uses hourly voice-traffic key performance indicators (KPIs) together with calendar covariates, and jointly predicts both KPIs over horizons 1:12 for each of ten cells. All real-data splits are chronological. Section IV-C gives the context lengths, split sizes, anchor strides, and rollout details.

TABLE I S YNTHETIC TASK PROTOCOLS . D EPLOYMENT TRAIN / TEST LENGTHS ARE COUNTED AFTER WASHOUT. F ORECAST HORIZONS ARE GIVEN IN TASK TIME STEPS . Task

Washout Deploy train

NARMA10

100

800

NARMA20

100

1200

NARMA30

100

1200

NARMA50

100

1600

MC

100

2000

Inubushi Lorenz63

100 200

1000 2000

SF-MG30

200

2000

MG84

200

1000

Lorenz96

200

2000

Test Input, target, and task parameter 800 ut ∼ U (0, 0.5); NARMA target, order k = 10 1200 ut ∼ U (0, 0.5); NARMA target, order k = 20 1200 ut ∼ U (0, 0.5); NARMA target, order k = 30 1600 ut ∼ U (0, 0.5); NARMA target, order k = 50 2000 ut ∼ U (−1, 1); delayed targets at lags 1:200 1000 ut ∼ U (−1, 1); yt = sin(2ut−5 ) 2000 scalar Lorenz63 forecast, (σ, ρ, β) = (10, 28, 8/3), horizon h = 25 2000 Mackey–Glass forecast, τ = 30, horizon h = 10 1000 Mackey–Glass forecast, τ = 17, horizon h = 84 2000 5D Lorenz96 forecast, forcing F = 8, horizon h = 25

b) Zero-rollout and exhaustive selectors: FP is the proposed zero-rollout selector (Algorithm 1). It ranks candidate operating points by deterministic kernel NMSE on labelled task-specific pilot data, without instantiating a finite reservoir. In the synthetic suite, three independently generated pilot sequences are used and each candidate is assigned its largest validation NMSE across the three pilots; this conservative aggregation penalizes pilot-sensitive configurations. The three pilot sequences play, for the deterministic FP selector, the role that reservoir seeds play for the simulation-based selectors: FP is replicated over independently sampled task pilots, whereas the empirical selectors are replicated over finite-reservoir realizations. Memoryn is a task-agnostic empirical proxy. At each grid point and selection width n, a ridge readout reconstructs D = 2n delayed samples of a generic scalar input ut ∼ U (−1, 1) from the reservoir features. Its score follows Jaeger’s memory-capacity criterion [28]: P   D (ut−d − ût−d )2 1 X MemScore(θ) = max 0, 1 − Pt , 2 D t (ut−d − ūd ) d=1 (32) where ūd is the sample mean of the lag-d target. Scores are averaged over three reservoir seeds. Directn evaluates the downstream validation task at width n and ranks candidates by their mean score over three reservoir seeds. ESN Directn uses the same task-direct protocol with the standard nonlinear ESN recurrence xt+1 = (1 − α)xt + α tanh(Wr xt + Win ut ). c) Budgeted selectors: We additionally compare three selectors under a matched retained candidate budget K. FP-K first ranks the admissible grid with FP at zero rollout cost and then applies Direct500 only to the top K candidates; we write FP-5, FP-10, etc. for specific budgets. Random-K samples K admissible candidates uniformly without replacement and returns the candidate with the best Direct500 validation score. TPE-K, a sequential Bayesian-optimization baseline, uses the outcomes of previous trials to direct later evaluations toward promising regions of the search space. Specifically, it applies the tree-structured Parzen estimator (TPE) [6] to the discrete

7

(σr , σin , α) search axes for K task-direct trials and returns the best observed candidate. Each task-direct trial uses three selection seeds, so one run of any budgeted selector spends 3K finite-reservoir rollouts; FP’s deterministic pre-ranking adds none. Random-K and TPE-K are repeated for 20 independent selection draws on the synthetic suite and 10 draws on ETT and Telco, reflecting the higher cost of the forecasting evaluations. Each individual run still receives the same matched budget of K candidates; repetitions are used only to estimate stochastic-selector variability. d) Regularization protocol: All ridge regressions in this paper (FP kernel selection, empirical selection, and deployment) use a single canonical regularization grid λ ∈ {10−12 , 10−11 , . . . , 102 } expressed on the kernel scale. Finitereservoir readouts are fitted on width-normalized features √ zt / n, so that the empirical Gram matrix n1 ZZ ⊤ and the deterministic FP kernel live on the same scale and a given λ corresponds to the same effective regularization at every stage and every reservoir width. The selection rule for λ differs by regime. On the synthetic benchmarks, where the data-generating processes are stationary and fresh i.i.d.-driven sequences can be generated, λ is selected by the holdout-ridge protocol on this grid. On the real-data ETT and Telco benchmarks, λ is selected by generalized cross-validation (RidgeCV) on the corresponding training block over the same grid, both inside finite-reservoir selection rollouts and at deployment. This uses the full training block for regularization selection, avoiding an additional dependence on a single chronological validation window in the nonstationary real-data setting. e) FP kernel implementation: All finite-reservoir deployments use coordinate-wise tanh features. For FP selection, we use the exact tanh kernel, evaluated by Gauss–Hermite quadrature, on all real-data benchmarks (ETT and Telco). For the synthetic benchmark suite and the synthetic timing sweeps, we use the closed-form erf surrogate from Remark 1. This surrogate is used only as a computational accelerator in the controlled synthetic setting; Appendix F-A verifies that matched erf-surrogate and exact-tanh selections have negligible selection-level impact on the checked synthetic settings. f) Synthetic benchmark protocol: The synthetic Cartesian grid contains 6 480 raw points. Applying the fixed stability guard (31) leaves 5 220 admissible candidates. Full-grid empirical selection at nselect = 500 therefore uses 5220×3 = 15 660 finite-reservoir rollouts per task, or 156 600 across the tentask suite. FP ranks the same admissible grid with zero finitereservoir rollouts. Its pilot sequences are generated after the task-specific warmup and have length Tpilot = 500, split into 333 training and 167 validation time steps; the deterministic selection kernel uses context L = 50. The washout lengths in Table I apply only to finite-reservoir deployment, since FP selection does not instantiate or roll out a reservoir. After selection, θ⋆ is fixed and evaluated with new reservoir seeds. Intermediate-width trends use 3 deployment seeds at n ∈ {1 000, 3 000, 5 000, 10 000}; the n = 20 000 synthetic comparison uses 10 deployment seeds. In the posthoc synthetic budget sweep, selected configurations use 3 deployment seeds.

g) Real-data benchmark protocol: All forecasting experiments use physical horizons of 1, 6, and 12 hours. A forecast anchor is a time index at which the preceding input history is used to predict the future target vector, and hence defines one forecasting example. Because nearby anchors in the long ETT series produce strongly overlapping examples, we retain every sth valid anchor, using the same anchor indices for all selectors. For each anchor, the finite-reservoir state is obtained by a local rollout over the input history represented by the FP kernel. In the FP construction, L is the maximum retained lag, so lags 0, . . . , L are included and each state uses L + 1 input samples. Thus, L = 96 corresponds to 97 hourly samples spanning a maximum lag of 96 hours, whereas L = 384 corresponds to 385 quarter-hourly samples over the same maximum-lag span. The FP kernel is evaluated only between the selected forecast anchors. For ETT, the input contains oil temperature (OT) and calendar covariates, and the target is future OT. ETTh is sampled hourly and uses L = 96, horizons h ∈ {1, 6, 12}, and anchor stride 43. ETTm is sampled every 15 minutes and uses L = 384, horizons h ∈ {4, 24, 48}, and anchor stride 173. Hence both variants use a 96-hour maximum retained lag and the same physical forecast horizons. ETTh uses 8 640/2 880/2 880 train/validation/test rows, while ETTm uses 34 560/11 520/11 520 rows. For the cellular-traffic task, the input combines data and voice KPIs with sine–cosine calendar covariates, and the target is the joint forecast of both KPIs over horizons 1:12. Each cell uses 841 training, 360 validation, and 932 test rows, with L = 96. Because the Telco series are shorter than ETT, every valid anchor is retained. In all forecasting experiments, one readout jointly predicts all horizons from a single reservoir state. Accordingly, one configuration θ⋆ is selected per ETT dataset or per Telco cell using the validation score aggregated over all horizons. Every selector, including full-grid Direct500 , optimizes this same aggregate objective. Selecting horizon-specific configurations is also possible by restricting the targets to one horizon, but is not considered here. The ETT and Telco Cartesian grids each contain 245 raw candidates, of which 154 satisfy the stability guard. All selectors operate on this same admissible set, so exhaustive empirical selection requires 154 × 3 = 462 finite-reservoir rollouts per ETT dataset and per Telco cell. For FP selection, the chronological training and validation blocks define the pilot split: the deterministic kernel model is fitted on the training anchors and evaluated on the validation anchors. Since each real-data series is fixed, a single pilot split is used rather than averaging over independently generated pilot sequences. All real-data models are deployed at width n = 20 000 with three fresh reservoir seeds. For ETT, reported deviations are computed over deployment seeds for deterministic FP-K and over independent selection draws for Random-K and TPE-K, after averaging deployment seeds within each draw. For Telco, results are first averaged over selection draws within each

8

TABLE II AVERAGE DEPLOYMENT SCORE ACROSS TEN SYNTHETIC BENCHMARKS UNDER THE HOLDOUT- RIDGE PROTOCOL . S CORES ARE 1 − NRMSE; HIGHER IS BETTER . M EAN ± STANDARD DEVIATION IS ACROSS TASKS . B EST AND SECOND - BEST RESULTS ARE BOLD AND UNDERLINED , RESPECTIVELY.

Selection rollouts

n = 1 000

n = 3 000

n = 5 000

n = 10 000

n = 20 000

0 0.692±0.171 0.744±0.159 0.756±0.158 0.768±0.158 0.772±0.158 156,600 0.619±0.194 0.668±0.193 0.683±0.197 0.698±0.188 0.703±0.186 156,600 0.720±0.178 0.749±0.164 0.762±0.158 0.769±0.157 0.774±0.156 156,600 0.708±0.203 0.725±0.195 0.727±0.188 0.731±0.188 0.734±0.186

cell and then summarized by the mean and standard deviation across the 10 cells. h) Metrics and statistical comparisons: Performance is reported as 1 − NRMSE; higher is better. The source of averaging and variability (tasks, deployment seeds, selection draws, or cells) is specified in the corresponding experimental protocol. Paired selector comparisons use Wilcoxon signedrank tests over tasks or cells, with Holm correction across the reported pairwise comparisons. Stochastic budget repetitions are used to characterise selection variability, not as additional independent tasks. B. Synthetic benchmark results The synthetic suite evaluates three complementary uses of the FP selector: zero-rollout selection when no simulation budget is available, width-portable deployment as the reservoir grows, and FP-based pre-screening when a small rollout budget can be spent. The paragraphs below address these in turn. a) Zero-rollout selection: Direct500 attains the highest mean deployment score across the ten-task suite, reaching 0.774 at n = 20 000 (Table II). This is expected: it directly evaluates finite-reservoir downstream validation performance during selection. FP reaches 0.772 while using zero finitereservoir search rollouts. The task-level comparison in Figure 1 shows that FP is best or tied-best on MC, NARMA30, NARMA50, and Lorenz96, and remains close to the best selector on several other tasks. The largest gaps occur on Lorenz63 and Inubushi, where finite-reservoir or ESN direct selection offers a clearer advantage. FP also exceeds ESN Direct500 on average at n = 20 000 (0.772 versus 0.734), showing that the analytically tractable linear-recurrence setting remains competitive in this suite. Full task-level scores are in Appendix E. b) Average-rank comparison: Figure 2 gives the rankbased comparison at n = 20 000. Direct500 has the best mean rank (1.65), followed by FP (2.05), ESN Direct500 (2.90), and Memory500 (3.40). After Holm correction at α = 0.05, FP and Direct500 are significantly better than Memory500 , while the pairwise differences among FP, Direct500 , and ESN Direct500 are not significant. The supported conclusion is that FP remains competitive with the task-direct simulation-based selectors while eliminating finite-reservoir search rollouts. c) Width portability and selection cost: The predictive comparison above is obtained under sharply different selection costs. Figure 3 reports measured full-grid selection time at several selection widths. FP does not take nselect as an algorithmic

Score (1 − NRMSE)

FP Memory500 Direct500 ESN Direct500

1.00

Score (1 − NRMSE)

Method

1.00

FP

Memory500

MC

NARMA10

NARMA20

NARMA30

NARMA50

Inubushi

Lorenz63

SF-MG30

MG84

Lorenz96

Direct500

ESN Direct500

0.75 0.50 0.25

0.75 0.50 0.25

Fig. 1. Deployment score at n = 20 000 across the ten-task synthetic suite under the holdout-ridge protocol. Scores are 1 − NRMSE; higher is better. Error bars show the standard deviation over deployment seeds.

4

Memory500 (3.40) ESN Direct500 (2.90)

3

2

1

Direct500 (1.65) FP (2.05)

Fig. 2. Average ranks across the ten synthetic tasks at n = 20 000. Lower is better. Horizontal cliques indicate no significant pairwise difference under Wilcoxon signed-rank tests with Holm correction (α = 0.05).

input: its cost depends on the pilot length, context length, candidate grid, and ridge grid, but not on reservoir width. We therefore summarize repeated FP timings by a horizontal mean line; the variation among the individual measurements reflects run-to-run timing variability rather than width dependence. The empirical selectors, by contrast, instantiate, roll out, and evaluate finite reservoirs at every candidate, so their measured cost increases with nselect ; the connecting lines are only guides to the eye. At deployment, the same FP operating point is reused without rerunning the selector. Its observed mean score increases from 0.692 at n = 1000 to 0.772 at n = 20 000, while the gap to Direct500 narrows from 0.028 to 0.002 (Table II). This narrowing is consistent with the large-width construction: as deployment width increases, the empirical feature geometry approaches the deterministic kernel targeted by FP. The result does not imply monotone improvement at every finite width. The task-level trends in Appendix E are consistent with this interpretation: FP is stable or improving over the reported widths on every task, whereas selectors tuned at nselect = 500

9

TABLE III B UDGETED SYNTHETIC COMPARISON AT n = 20 000. ROWS INDICATE THE RETAINED CANDIDATE BUDGET K; SCORES ARE DEPLOYMENT 1 − NRMSE ( HIGHER IS BETTER ), AVERAGED ACROSS THE TEN TASKS . P ERCENTAGES ARE RELATIVE TO EXHAUSTIVE D IRECT 500 SELECTION .

K 25 50 100 250

Ten-task rollouts

FP-K

can exhibit non-monotone trends on individual tasks. d) Budgeted screening: The zero-rollout selector can also serve as a pre-screen when a small simulation budget is available. FP first ranks the full grid without finite-reservoir rollouts; Direct500 is then evaluated only on the top K candidates. Table III compares this FP-K strategy with RandomK and TPE-K, using the same retained candidate budget but reporting Random/TPE over 20 draws to estimate baseline variability. FP-K gives the strongest mean performance at every reported budget. The advantage is largest at small K and shrinks as K grows, indicating convergence of the budgeted selectors. After Holm correction across the eight FP-vs-baseline budget comparisons, the paired differences are directional rather than significant. The main advantage of FP-K is therefore budget efficiency and determinism: at K = 250, FP-K matches the full-grid Direct500 reference within sampling variation (0.777 versus 0.774) using 750 finite-reservoir selection rollouts per task (4.8% of the full grid). Random-K and TPE-K use the same per-run candidate budget, but their reported means are estimated over 20 stochastic draws and reach 0.750 and 0.769, respectively. The screened protocol therefore provides an intermediate operating mode between zero-rollout selection and exhaustive task-direct search, and motivates the low-budget Telco experiments (FP-5, FP-10) below. e) Selection-length sensitivity: An auxiliary coarse-grid sweep over the selection-sequence length (Appendix F) shows that FP is close to its highest peak already at Tselect = 250 and reaches 0.771 by Tselect = 500, whereas direct selection continues to benefit from longer sequences as additional taskspecific validation data provide a more reliable ranking of the candidate reservoirs. Memory500 remains constant because its task-agnostic memory proxy selects the same candidate at every tested selection length. FP is therefore efficient not only in finite-reservoir rollouts but also in the length of the labelled pilot sequence required for selection.

TPE-K

750 (0.5%) 0.774±0.159 0.669±0.234 0.660±0.244 1,500 (1.0%) 0.779±0.162 0.701±0.232 0.716±0.214 3,000 (1.9%) 0.775±0.168 0.732±0.219 0.742±0.209 7,500 (4.8%) 0.777±0.166 0.750±0.194 0.769±0.167

Full Direct500 reference

Fig. 3. Measured aggregate selection CPU time over the ten synthetic tasks as a function of selection width nselect . All finite-reservoir selector points are measured; connecting lines are guides to the eye and do not represent fitted complexity laws. FP does not take nselect as an input, so its repeated measured runtimes are summarized by a horizontal mean line, with individual markers and a shaded ±1 standard-deviation band showing run-to-run timing variability. FP uses the closed-form erf surrogate employed for synthetic selection; see Sec. III-C. Lower is better.

Random-K

0.774±0.156

C. Real-data benchmark results We complement the synthetic suite with real-data forecasting experiments derived from five datasets: the four public ETT-small datasets (ETTh1, ETTh2, ETTm1, and ETTm2) and an operational cellular-traffic dataset for ten randomly selected network cells. The corresponding protocol details are given in the real-data benchmark protocol of Section IV. The ETT datasets provide reproducible public benchmarks, while the Telco task tests the selector in an operational multivariate setting where finite-reservoir validation rollouts are costly. We therefore focus on two questions: whether zero-rollout FP selection remains competitive with exhaustive task-direct search, and whether FP-based screening allocates small rollout budgets more effectively than Random-K and TPE-K. a) Zero-rollout selection: Table IV reports datasetaverage scores for all selectors; horizon-level results are in Appendix G. On ETTm1, zero-rollout FP is the strongest of the linear-recurrence selectors (0.505 versus 0.499 for fullgrid Direct500 ), while ESN Direct500 is slightly higher (0.508). At the longest ETTm horizons, the task-agnostic Memory500 proxy is the strongest selector (Appendix G). This reflects a short-versus-long-horizon trade-off under aggregate selection. On ETTm2, for example, Memory500 improves the 6- and 12-hour scores to 0.474 and 0.501, compared with 0.456 and 0.442 for FP, but reduces the 1-hour score from 0.862 to 0.762. On the ETT dataset averages, Memory500 remains below the task-informed selectors throughout (Table IV). These results also indicate that the ETT candidate landscape is relatively flat: several operating points have nearly identical deployment scores, so stochastic search can catch up as the budget grows. Telco is different. On the three representative horizons shown in Table IV, zero-rollout FP trails full-grid Direct500 (0.572 versus 0.635), and the task-agnostic Memory500 proxy is also higher (0.619). This shows that zero-rollout ranking alone is not uniformly sufficient on real data. However, the budgeted aggregate comparison in Table V shows that only a very small task-direct budget is needed to close the gap: FP-5 reaches 0.622 ± 0.112 against the full-grid Direct500 reference 0.621 ± 0.113, using 15 rather than 462 selection rollouts per cell. b) Budgeted screening: As on the synthetic suite, FP-K ranks the grid deterministically at zero rollout cost and spends the task-direct budget only on the retained candidates. Table V compares FP-K with Random-K and TPE-K at matched candidate budgets: each individual selector evaluates K candidates

10

TABLE IV T HREE - HORIZON REAL - DATA FORECASTING SUMMARY AT n = 20 000. E ACH ENTRY IS THE MEAN±STANDARD DEVIATION OVER THE THREE REPORTED PHYSICAL HORIZONS (1, 6, AND 12 H ). H ORIZON - LEVEL SCORES , WHOSE DEVIATIONS ARE INSTEAD OVER DEPLOYMENT SEEDS (ETT) OR CELLS (T ELCO ), ARE GIVEN IN S EC . S-VII OF THE SM. T HE T ELCO AGGREGATE OVER HORIZONS 1:12 IS REPORTED IN TABLE V. S CORES ARE 1 − NRMSE; HIGHER IS BETTER .

Dataset

FP

Memory500

Direct500

ESN Direct500

ETTh1 0.556±0.16 0.522±0.11 0.567±0.17 ETTh2 0.664±0.15 0.525±0.05 0.664±0.15 ETTm1 0.505±0.19 0.498±0.18 0.499±0.19 ETTm2 0.587±0.24 0.579±0.16 0.587±0.24 Telco 0.572±0.03 0.619±0.05 0.635±0.04

0.564±0.18 0.662±0.15 0.508±0.20 0.587±0.24 0.607±0.03

using three reservoir seeds, hence 3K finite-reservoir rollouts. Random-K and TPE-K are repeated for ten independent draws to estimate their variability; these repetitions use 30K reported rollouts in total but do not increase the candidate budget available to any individual selector. Their standard deviations therefore quantify selection instability: a single stochastic run has the same 3K budget but can be substantially below the reported mean. On the public ETT benchmarks, FP-K recovers the fullgrid Direct500 operating point already at K = 5 on ETTh1, ETTh2, and ETTm2, and selects the same operating point at every reported budget. TPE-K attains its strongest ETTh1 test score at K = 10, while the three screening methods converge at larger budgets. On ETTh2, Random-K and TPE-K improve with K but remain below FP-K even at K = 50; on ETTm2 they trail at small budgets and approach FP-K only at larger K. The operational Telco task is more discriminative and provides the strongest evidence for pilot-informed screening. Full-grid Direct500 reaches 0.621 ± 0.113 across cells, at 154 × 3 = 462 rollouts per cell. FP-5 already matches or slightly exceeds this reference (0.622±0.112) with 15 rollouts, compared with 0.548 for Random-K and 0.572 for TPE-K at the same budget, and FP-K remains between 0.621 and 0.623 at every reported budget, compared with 0.584 and 0.593 for the baselines at K = 10. FP-K beats both baselines on all ten cells at K = 5, K = 10, and K = 25; after Holm correction over the eight FPvs-baseline budget comparisons, the adjusted values satisfy pHolm ≤ 0.0078 in all three cases. FP-50 recovers the full-grid Direct500 operating point with 150 rather than 462 rollouts, and by K = 50 the methods have largely converged and the corrected differences are no longer significant. Taken together, the forecasting results show that FP ranking is most valuable when validation rollouts are scarce. It either recovers the full-grid operating point immediately, as on several ETT datasets, or allocates a small task-direct budget to substantially better candidates than Random-K and TPE-K, as on the operational Telco task. V. D ISCUSSION a) Width-independent tuning of reservoir dynamics: The large-width kernel relocates selection cost from the reservoir to

the pilot sequence. Empirical selectors instantiate a reservoir for every grid point; FP instantiates none. This has two concrete consequences. First, the selection cost is independent of reservoir width: the same θ⋆ computed at pilot length Tpilot = 500 can be deployed at any n without rerunning selection. In our experiments, its mean score improves over the reported widths, consistently with convergence toward the large-width kernel. Second, cost scales with pilot length rather than with reservoir width. The selection-length study (Appendix F) shows that FP plateaus at Tselect = 500, while task-direct selection continues to improve; FP is therefore dataefficient as well as simulation-efficient. b) Deterministic screening versus stochastic search: Random-K and TPE-K are important baselines because they represent the finite-rollout strategies a practitioner would normally try instead of full-grid selection. Their behavior also clarifies the role of FP. FP is not a replacement for high-budget sequential search; rather, it provides a deterministic ranking of the full candidate grid before any reservoir is instantiated. The mechanism behind the low-budget advantage is that the FP ranking is task-informed from the outset: the pilot kernel already encodes the temporal structure of the downstream task, so the first K candidates evaluated are concentrated in the relevant region of the grid. Random-K carries no task information, and TPE-K must spend its earliest trials learning the response surface before it can exploit it. The budgeted results show that this ranking is most valuable at small K. On the synthetic suite, FP-K gives the strongest mean score at every reported budget and matches the full-grid Direct500 reference within sampling variation using 4.8% of the full-grid rollout budget. On Telco, exact-tanh FP-K is already at the full-grid Direct500 level with K = 5 and remains stable across the reported budgets, while Random-K and TPE-K require larger budgets to approach the same operating point. The lowbudget differences against both stochastic baselines are significant for K ≤ 25 after Holm correction (pHolm ≤ 0.0078). On flatter ETT landscapes, the stochastic baselines catch up as K grows, but their reported means are expectations over stochastic selector runs; a single run receives the same 3K candidate budget and can vary substantially, as reflected by the standard deviations. Thus, the practical distinction is not only cost, but also determinism and stability of the selected operating point. c) Recurrent topology as a kernel design variable: The theory requires two properties of the recurrent matrix: asymptotic operator-norm control and diagonal self-averaging of the entries of Akn (Aℓn )⊤ (Lemma 1). Ginibre matrices satisfy both with probability tending to one. The same framework extends to cyclic-shift and Haar-orthogonal recurrences, and more generally to ensembles with deterministic crosslag propagation coefficients; the corresponding statements and coefficient catalogue are provided in Appendices C–D. In this sense, the recurrent matrix is not merely an implementation detail: once its limiting memory coefficients are specified, it becomes part of the analytic kernel used for selection. d) Scope and extensions: The exact theory requires linear recurrent dynamics. For a nonlinear ESN, the state is no longer a linear combination of propagated inputs, and the large-

11

TABLE V B UDGETED REAL - DATA FORECASTING COMPARISON AT n = 20 000. E ACH SELECTOR EVALUATES K CANDIDATES USING 3K FINITE - RESERVOIR ROLLOUTS . † DENOTES FP-K SIGNIFICANTLY OUTPERFORMING BOTH R ANDOM -K AND TPE-K AFTER H OLM CORRECTION .

Dataset

Selector

K=5 15 rollouts

K = 10 30 rollouts

K = 25 75 rollouts

K = 50 150 rollouts

ETTh1

FP-K Random-K TPE-K FP-K Random-K TPE-K FP-K Random-K TPE-K FP-K Random-K TPE-K FP-K Random-K TPE-K

0.567±0.00035 0.561±0.008 0.561±0.013 0.665±0.00049 0.609±0.047 0.626±0.039 0.501±0.00079 0.493±0.013 0.493±0.014 0.587±0.00020 0.579±0.014 0.569±0.011 0.622±0.112† 0.548±0.170 0.572±0.144

0.567±0.00035 0.563±0.009 0.571±0.007 0.665±0.00049 0.645±0.015 0.635±0.032 0.499±0.00086 0.495±0.008 0.497±0.008 0.587±0.00020 0.581±0.012 0.571±0.008 0.622±0.110† 0.584±0.142 0.593±0.120

0.567±0.00035 0.568±0.002 0.567±0.005 0.665±0.00049 0.654±0.011 0.639±0.026 0.499±0.00086 0.498±0.009 0.503±0.006 0.587±0.00020 0.579±0.004 0.583±0.005 0.623±0.110† 0.598±0.136 0.602±0.127

0.567±0.00035 0.567±0.002 0.567±0.005 0.665±0.00049 0.661±0.003 0.660±0.007 0.499±0.00086 0.499±0.007 0.501±0.004 0.587±0.00020 0.584±0.002 0.584±0.003 0.621±0.113 0.617±0.116 0.612±0.114

ETTh2

ETTm1

ETTm2

Telco

width covariance involves time-varying activation factors that depend on the input trajectory. A frozen-diagonal or meanfield approximation would replace the exact trace moments with approximate ones, at the cost of the exactness guarantee. Two further limitations delimit the intended use of the selector. First, the pilot sequence must be drawn from the same task distribution as deployment: FP ranks operating points for the task the pilot describes, and its behavior under pilot–deployment distribution shift has not been studied here. Second, FP ranks a fixed discrete candidate grid and does not refine candidates continuously; a natural hybrid uses FP to screen the grid and a sequential method to refine locally within the retained region. Extending the framework to nonlinear recurrences, deep or stacked reservoir layers, and structured state-space models are natural directions for future work. VI. C ONCLUSION We introduced a deterministic zero-rollout hyperparameter selector for reservoir-based temporal sequence learning. The selector uses a free-probability large-width kernel: the Gram matrix of nonlinear temporal features converges to a deterministic limit computable from a short labelled pilot sequence and the cross-lag propagation coefficients of the reservoir ensemble. This moves hyperparameter search offline: no reservoir is instantiated during selection, the cost is independent of deployment width, and the selected configuration transfers across widths without rerunning the search. Across ten temporal benchmarks covering nonlinear sequence processing, memory, and chaotic forecasting, Holmcorrected pairwise tests find no significant difference between FP and the task-direct simulation-based selectors requiring 156 600 finite-reservoir selection rollouts, while FP is significantly better than the task-agnostic Memory baseline. In the budgeted synthetic sweep, FP-K gives the strongest mean score at each reported budget and matches the full-grid Direct500 reference within sampling variation using only 4.8% of the full-grid rollout budget. On public ETT data, FP-K recovers

Full Direct500 0.567±0.00035

0.665±0.00049

0.499±0.00086

0.587±0.00020

0.621±0.113

the full-grid Direct500 operating point on three of four datasets already at K = 5. On the ten-cell Telco forecasting task, FP-5 slightly exceeds full-grid Direct500 (0.622 versus 0.621) using 15 rather than 462 selection rollouts per cell, and FP-K is significantly better than Random-K/TPE-K for K ≤ 25 after Holm correction (pHolm ≤ 0.0078). The results support free-probability kernels as a practical tool for pilot-informed reservoir model selection: they are most useful when finite-reservoir validation rollouts are scarce, while remaining compatible with higher-budget direct or sequential search. Future work will extend the construction to nonlinear recurrences and broader classes of structured sequence models. ACKNOWLEDGMENTS This work was supported by the Norwegian Research Council SFI NorwAI (309834); by EMERGE, a project funded by the European Innovation Council (prj. code 101070918); and by NEURONE, a project funded by the European Union – Next Generation EU, M4C1 CUP I53D23003600006, under program PRIN 2022 (prj. code 20229JRTZA). A PPENDIX A F IXED -C ONTEXT K ERNEL A. Proof roadmap and notation The proof follows four steps. First, the mixed propagation traces converge to deterministic coefficients τk,ℓ . Second, the corresponding diagonal entries converge to the same coefficients on average over reservoir coordinates. Third, conditional on the recurrent matrix, each coordinate pair is Gaussian; the diagonal result therefore identifies the limiting nonlinear feature kernel. Finally, concentration with respect to the Gaussian input matrix removes the remaining finite-width fluctuations. We use the same notation as in the main text. For each width n, let σr , An = (1 − α)In + αWr(n) , Wr(n) = √ G(n) n r

12

(n)

where the entries of Gr are independent standard Gaussian variables. Independently, σin (n) (n) Win = √ Gin , din (n)

with Gin again standard Gaussian. For deterministic u, v ∈ Rdin , h i σ2 (n) (n) (33) EWin Win uv ⊤ (Win )⊤ = in (u⊤ v)In . din For a fixed maximum retained lag L, define (L)

xt

:= α

L X

(n)

Akn Win ut−1−k .

(34)

k=0

Thus, lags 0, . . . , L are retained, for up to L + 1 input time points. For fixed k, ℓ ≥ 0, set min{k,ℓ} 

τk,ℓ :=

X j=0

k j

  ℓ (1 − α)k+ℓ−2j (ασr )2j . j

(35)

The deterministic linear-state covariance and its associated 2× 2 matrix are σ2 (L) Qθ (t, s) := α2 in din

L X

τk,ℓ u⊤ t−1−k us−1−ℓ ,

(L) Σθ (t, s) :=

(L)

Qθ (t, t) (L) Qθ (t, s)

! (L) Qθ (t, s) . (L) Qθ (s, s)

(37)

For a coordinate-wise feature map ψ, define (L)

Kψ,n,θ (t, s) :=

1 (L) ψ(xt )⊤ ψ(x(L) s ), n

(38)

and (L)

Kψ,θ (t, s) := E[ψ(gt )ψ(gs )] .

where mn (i, j) := n−1 Tr(Xni (Xn⊤ )j ). Normalized real Ginibre matrices converge in expected ∗moments to a standard circular element c [34, Sec. 11.6.3]. The free Wick rule gives
φ ci (c∗ )j = δij : a non-crossing pairing of the ordered word ci (c∗ )j exists only when i = j, and is then unique. Hence E[mn (i, j)] → δij . It remains to control fluctuations. Differentiating mn (i, j) with respect to the Gaussian entries and summing the resulting position terms yields  (i + j)2 max 1, ∥Xn ∥2(i+j−1) . ∥∇mn (i, j)∥2 ≤ 2 n The Gaussian Poincaré inequality and (40) therefore give Var(mn (i, j)) = O(n−2 ). Thus mn (i, j) → δij in L2 . Substitution into the finite sum (42) gives (41) with the coefficients (35).

(36)

k,ℓ=0

and

(n) √ Proof. For the proof, write Xn := Gr / n. Since In commutes with X
n , with the binomial theorem we can write n−1 Tr Akn (Aℓn )⊤ as k X ℓ    X k ℓ (1 − α)k+ℓ−i−j (ασr )i+j mn (i, j), (42) i j i=0 j=0

(39)

(L)

Here (gt , gs ) ∼ N (0, Σθ (t, s)). Throughout, the input sequence is deterministic and bounded: U := supt ∥ut ∥ < ∞. The symbol C denotes a finite constant that may change from line to line, may depend on fixed model parameters, L, U , and ψ, but never on n. We use repeatedly the standard normalized-Ginibre bound # " (n) m Gr √ sup E <∞ for every fixed m. (40) n n

C. Diagonal self-averaging Lemma 2 (Averaged diagonal self-averaging). For every fixed L < ∞, n 1 X  k ℓ ⊤ 2 P → 0. (43) max An (An ) ii − τk,ℓ − 0≤k,ℓ≤L n i=1 Proof. Fix k, ℓ and set Bn := Akn (Aℓn )⊤ . Simultaneous (n) permutation of the rows and columns of Gr preserves its law. Hence the diagonal entries of Bn are identically distributed and   1 E[(Bn )ii ] = E Tr Bn −→ τk,ℓ n by Proposition 2. To control one diagonal entry, expand Bn into the finitely (n) √ many matrix words in Xn := Gr / n and Xn⊤ . For a word M = F1 · · · Fr , differentiation of (M )11 at each occurrence of Xn or Xn⊤ gives X

∂(M )11

a,b

(n) ∂(Gr )ab

2

≤

 r2 max 1, ∥Xn ∥2(r−1) . n

Indeed, the expected norm is at most 2 + o(1) and Gaussian concentration gives a sub-Gaussian upper tail [33]. Consequently, every fixed power of An has uniformly bounded norm moments.

After summing over the finitely many words, Gaussian Poincaré and (40) yield Var((Bn )11 ) = O(n−1 ). Therefore # " n   1X 2 E |(Bn )ii − τk,ℓ | = E |(Bn )11 − τk,ℓ |2 −→ 0. n i=1

B. Mixed propagation moments

Markov’s inequality and a finite union bound over 0 ≤ k, ℓ ≤ L prove (43).

Proposition 2 (Ginibre recurrent memory coefficients). For every fixed k, ℓ ≥ 0,
L1 1 Tr Akn (Aℓn )⊤ −−→ τk,ℓ . n

(41)

Only the diagonal propagation terms are needed to identify the marginal Gaussian law of each coordinate pair entering the kernel. Dependence between distinct reservoir coordinates is controlled separately by the concentration step below.

13

D. Nonlinear feature-kernel limit Theorem 2 (Deterministic nonlinear reservoir kernel). Let ψ ∈ Cb1 (R). For fixed L, t, s and a deterministic bounded input sequence,

A PPENDIX B C OMPLETE H ISTORY AND S ELECTION C ONSISTENCY

The convergence is with respect to the joint randomness of (n) (n) Wr and Win .

The fixed-context result is the object used directly by the selector. This section shows that, in the stable regime, increasing the maximum retained lag recovers the complete washedout history. The argument has three parts: a deterministic covariance-tail bound, a uniform finite-width power bound for Ginibre recurrence, and a final triangle-inequality passage to the nonlinear kernel.

Proof. Write Mψ := ∥ψ∥∞ and Lψ := ∥ψ ′ ∥∞ . For a positivesemidefinite 2 × 2 matrix Σ, let

A. Deterministic complete-history kernel

(L)

P

(L)

Kψ,n,θ (t, s) − → Kψ,θ (t, s)

h(Σ) := E[ψ(g1 )ψ(g2 )] ,

(n → ∞).

(44)

(g1 , g2 ) ∼ N (0, Σ).

Coupling two Gaussian vectors through a common standard normal vector and using ∥Σ1/2 − (Σ′ )1/2 ∥ ≤ ∥Σ − Σ′ ∥1/2 [35, Thm. X.1.1] gives √ (45) |h(Σ) − h(Σ′ )| ≤ 2 2 Mψ Lψ ∥Σ − Σ′ ∥1/2 .

Proposition 3 (Complete-history covariance and context limit). Assume σr < 1 and U = supt ∥ut ∥ < ∞. Define σ2 X (∞) τk,ℓ u⊤ (50) Qθ (t, s) := α2 in t−1−k us−1−ℓ . din k,ℓ≥0

The series converges absolutely. With rθ := 1 − α + ασr2 < 1,

(L) (L) Conditional on An , the pair (xt,i , xs,i ) is centered Gaus-

sian. Its covariance entries are (L)

cn,i (r, q) = α

2 2 σin

din

L X

(∞)

Qθ

 k ℓ ⊤ ⊤ An (An ) ii ur−1−k uq−1−ℓ . (46)

(L)

(t, s) − Qθ (t, s) ≤

(∞)

(∞)

If Kψ,θ is obtained from Qθ

k,ℓ=0

(L) Let Cn,i (t, s) be the corresponding 2 × 2 covariance matrix.

By Lemma 2,

(∞)

(L)

(L+1)/2

.

(52)

(L)

(∞)

(L → ∞).

(53)

Kψ,θ (t, s) − Kψ,θ (t, s) ≤ Cθ rθ In particular,

Proof. The coefficients τk,ℓ are nonnegative and have generating function X F (x, y) := τk,ℓ xk y ℓ (47)

It remains to control fluctuations around the conditional (n) mean. Stack the entries of Gin into a standard Gaussian vector (L) (L) g and write xt = Mt g, xs = Ms g. Since

k,ℓ≥0

=

1 . (54) (1 − (1 − α)x)(1 − (1 − α)y) − α2 σr2 xy

At x = y = 1, X

L X

ασin U ∥Mt ∥ ≤ √ ∥Akn ∥, din k=0 the conditional Gaussian Poincaré inequality gives !2 L   C X (L) k Var Kψ,n,θ (t, s) | An ≤ ∥An ∥ . n

τk,ℓ =

k,ℓ≥0

(48)

k=0

1 , α2 (1 − σr2 )

which proves absolute convergence for bounded inputs. Moreover, X 1 1 F (x, 1) = , τk,ℓ = rθk . α(1 − rθ x) α ℓ≥0

2

The right-hand side is OP (n−1 ) by (40). Hence the conditional fluctuation vanishes in probability, and (47) proves (44). For reference, the same conditional-mean calculation before applying ψ, together with the Gaussian quadratic-form variance bound, also gives 1 (L) ⊤ (L) P (L) (x ) xs − → Qθ (t, s). n t

(51)

as in (39), then

Kψ,θ (t, s) −→ Kψ,θ (t, s)

n

2 P 1X (L) (L) Cn,i (t, s) − Σθ (t, s) − → 0. n i=1 F

Together with (45), this implies h i P (L) (L) E Kψ,n,θ (t, s) | An − Kψ,θ (t, s) − → 0.

2 2σin U2 rL+1 . din (1 − σr2 ) θ

(49)

On any fixed finite pilot set, a union bound gives entrywise and hence Frobenius-norm convergence of the associated kernel blocks.

The complement of {0, . . . , L} is contained in {k > L} ∪ {ℓ > L}. By symmetry, X 2 X k 2 τk,ℓ ≤ rθ = 2 rL+1 . α α (1 − σr2 ) θ max{k,ℓ}>L

k>L

2 Multiplying by α2 σin U 2 /din gives (51). The same bound (∞) (L) applies to each entry of Σθ −Σθ . The Gaussian continuity

estimate (45) then yields (52). Because the admissible candidate grid and pilot index sets are finite, the convergence in (53) is uniform over all candidates and pilot-kernel entries used by the selector.

14

Proposition 4 bounds this term geometrically in L, uniformly for all sufficiently large n. For fixed L, Theorem 2 gives

B. Finite-width Ginibre tail Proposition 4 (Uniform Ginibre power bound). Let σr An = (1 − α)In + α √ G(n) n r

(L)

(L)

and assume σr < 1. For any 1 < R < σr−1 , set qR := 1 − α + ασr R < 1. Almost surely, there are finite random constants CR and NR such that, for all n ≥ NR and m ≥ 0, m ∥Am n ∥ ≤ CR qR .

(55)

Consequently, for bounded input histories, (∞)

xt

:= α

∞ X

(n)

Akn Win ut−1−k

(56)

k=0

is well defined for all sufficiently large n, and αCR σin U 1 (∞) (L) ≤√ q L+1 lim sup √ xt − xt n din (1 − qR ) R n→∞

(57)

almost surely. (n) √ Proof. Write Xn := Gr / n for this proof. The normalized real Ginibre matrix converges strongly to a standard circular element c [36]. A standard consequence of strong convergence is uniform resolvent control on compact subsets of C\spec(c): by hermitization, the smallest singular value of zIn − Xn is eventually bounded away from zero uniformly on such compact sets. Since spec(c) is the closed unit disk [34], for every R > 1 there are almost surely finite MR , NR such that

sup sup

(zIn − Xn )−1 ≤ MR .

P

(L)

Kψ,n,θ (t, s) − → Kψ,θ (t, s),

(58)

n≥NR |z|=R

With f (z) = 1 − α + ασr z, we have An = f (Xn ). The holomorphic functional calculus and (58) give Z 1 ∥Am ∥ ≤ |f (z)|m (zIn − Xn )−1 |dz| n 2π |z|=R

(∞)

while Proposition 3 gives Kψ,θ (t, s) → Kψ,θ (t, s). The result follows from the triangle inequality by first choosing L large and then letting n → ∞. D. Selection consistency Corollary 2 (Consistency of FP selection). Assume that the deterministic kernel and the finite reservoir use the same feature map. Fix finite candidate and ridge grids, with all ridge parameters strictly positive. For any fixed maximum retained lag L, if the deterministic validation score has a unique minimizer, then finite-width reservoir selection returns the same minimizer with probability tending to one as n → ∞. In the stable regime, a unique complete-history optimum is preserved by all sufficiently large finite contexts, and selection based on the complete finite-reservoir trajectory converges to this optimum as n → ∞. Proof. Let T and V be the fixed pilot training and validation indices. For λ ≥ λmin > 0, the prediction map (KT T , KV T ) 7−→ KV T (KT T + λI)−1 YT is continuous and uniformly Lipschitz on the fixed finite blocks. This follows from ∥(KT T + λI)−1 ∥ ≤ λ−1 min and the resolvent identity. The same is therefore true of the validation NMSE, assuming its denominator is nonzero. Theorem 2 gives convergence of every fixed-context kernel entry. A union bound over the finite pilot, candidate, and ridge grids yields (L)

(L)

P

max NMSEval,n (θ, λ) − NMSEval (θ, λ) − → 0. θ,λ

C. Complete-history kernel

A unique minimizer on a finite grid has a positive gap to the second-best score. Whenever the uniform error is less than half this gap, the empirical and deterministic minimizers agree. This proves the fixed-context statement. Proposition 3 gives uniform convergence of the finitecontext deterministic scores to their complete-history counterparts over the same finite grids. Hence a unique completehistory minimizer is also the finite-context minimizer for all sufficiently large L. Finally, Theorem 3 gives uniform convergence of the empirical complete-history scores, and the same positive-gap argument proves the final claim.

Theorem 3 (Complete-history Ginibre kernel). Under the assumptions of Theorem 2 and Proposition 4, for every fixed t, s,

Near ties have a small score gap and may therefore require larger reservoir width or maximum lag before the ranking stabilizes.

m ≤ RMR qR .

Thus (55) holds with CR = RMR . Summing the geometric tail and using (n)

∥Win ∥ σin √ −→ √ n din

almost surely

proves (57).

1 P (∞) (∞) )− → Kψ,θ (t, s). ψ(xt )⊤ ψ(x(∞) s n Proof. By boundedness and Lipschitz continuity of ψ, (∞)

Kψ,n,θ (t, s) :=

(∞)

(59)

(L)

Kψ,n,θ (t, s) − Kψ,n,θ (t, s)  Mψ Lψ  (∞) (L) ∥xt − xt ∥ + ∥x(∞) − x(L) ≤ √ s s ∥ . n

(60)

A PPENDIX C H AAR -O RTHOGONAL R ECURRENT M ATRIX The Ginibre proof uses Gaussian concentration to establish diagonal self-averaging. For a Haar-orthogonal recurrent matrix, the same conclusion follows directly from orthogonal invariance.

15

Proposition 5 (Haar-orthogonal extension). Let On be Haar (n) distributed on the orthogonal group, independently of Win , and let An = (1 − α)In + ασr On . For every fixed L < ∞, Lemma 2 holds with the coefficients τk,ℓ in (35). Consequently, Theorem 2 holds with the same deterministic kernel. If σr < 1, the complete-history conclusion also holds. Proof. Fix k, ℓ and set Bn = Akn (Aℓn )⊤ . Orthogonal conjugation preserves the law of Bn . If v is independent and uniform on the unit sphere, then d

⊤

(Bn )11 = v Bn v. The spherical fourth-moment identity gives # " 2 2∥Bn ∥2 1 ⊤ . E v Bn v − Tr Bn Bn ≤ n n+2 Moreover, ∥Bn ∥ ≤ (1 − α + ασr )k+ℓ . The normalized traces of fixed Laurent polynomials in On converge to the corresponding Haar-unitary moments [37], which give the same τk,ℓ as in (35). Thus h i 2 E |(Bn )11 − τk,ℓ | −→ 0. The diagonal entries are exchangeable, so averaging over them, followed by Markov’s inequality and a finite union bound, proves the diagonal self-averaging statement. The fixedcontext kernel follows from the same conditional-Gaussian argument as Theorem 2. For σr < 1,

(L)

Here Qθ

is formed with the relevant coefficients ηk,ℓ . Since √ δerf := sup tanh(x) − erf( πx/2) < 0.0354, x

the corresponding limiting kernel entries differ by less than 2δerf < 0.071. The surrogate is used only for the synthetic selection and timing experiments; all ETT and Telco results use the exact tanh kernel evaluated by Gauss–Hermite quadrature. The empirical validation is reported in Sec. F-A. A. Ginibre, Haar-orthogonal, and cyclic-shift matrices For the real Ginibre matrix used in the experiments, σr Wr(n) = √ G(n) , n r Proposition 2 gives ηk,ℓ = τk,ℓ . Proposition 5 shows that (n) the same coefficients hold for Wr = σr On with On Haar orthogonal. They also hold for the deterministic cyclic-shift matrix Cn , defined by Cn ej = ej+1 (mod n) ,

Wr(n) = σr Cn .

Indeed, Cn⊤ = Cn−1 and, for fixed i, j and all sufficiently large n,
1 Tr Cni (Cn⊤ )j = δij . n The three structures therefore produce the same fixed-context coefficients. Their complete-history control differs: Ginibre uses Proposition 4, whereas the orthogonal matrices satisfy the direct geometric bound

m ∥Am n ∥ ≤ (1 − α + ασr ) ,

m ∥Am n ∥ ≤ (1 − α + ασr ) .

which gives the complete-history extension directly. B. Circulant matrices

A PPENDIX D OTHER R ECURRENT M ATRIX S TRUCTURES The kernel derivation depends on the recurrent matrix only through the limits of
1 Tr Akn (Aℓn )⊤ n and the corresponding diagonal self-averaging property. This section records these coefficients for several recurrent matrix structures. Once they are known, the covariance and nonlinear kernel are obtained from (36)–(39) by replacing τk,ℓ with the appropriate coefficient ηk,ℓ . No additional kernel notation is needed.

where Fn is the unitary discrete Fourier matrix. Assume that the empirical distribution of the uniformly bounded Fourier eigenvalues converges weakly to ν. Then Z circ ηk,ℓ = (1 − α + ασr z)k (1 − α + ασr z)ℓ dν(z). (62)

Remark 1 (Optional closed-form surrogate). For any of the real-valued recurrent matrices below, the choice √ ψerf (x) = erf( πx/2)

Then

yields a closed-form kernel. With (L)

Cn = Fn∗ diag(γn,0 , . . . , γn,n−1 )Fn ,

Wr(n) = σr Cn ,

C

If the matrices are generated by a limiting Fourier symbol g : [0, 2π] → C, define aθ (ω) = 1 − α + ασr g(ω). circ ηk,ℓ =

1 2π

Z 2π

ℓ

aθ (ω)k aθ (ω) dω.

0

Because Akn (Aℓn )⊤ is circulant, every diagonal entry equals its normalized trace; diagonal self-averaging is therefore automatic. Complete history follows whenever

(L)

Dθ (t) := 2 + πQθ (t, t), the kernel is   (L) πQ (t, s) 2 (L) θ . Kerf,θ (t, s) = arcsin q π (L) (L) Dθ (t)Dθ (s)

Let

sup |aθ (ω)| < 1. (61)

ω∈[0,2π]

The cyclic shift is the special case g(ω) = eiω .

16

TABLE VI R ECURRENT- MATRIX CATALOGUE . T HE COEFFICIENT COLUMN GIVES THE QUANTITIES SUBSTITUTED FOR τk,ℓ IN (36); THE NONLINEAR KERNEL IS THEN FORMED AS IN (39). Recurrent matrix

Coefficients

Complete-history control

Nonlinear-kernel status

Real Ginibre

ηk,ℓ = τk,ℓ

Resolvent bound

Haar orthogonal Cyclic shift Circulant Gaussian skew-symmetric Complex Ginibre / Haar unitary

ηk,ℓ = τk,ℓ ηk,ℓ = τk,ℓ for fixed lags and large n Fourier integral (62) Semicircle integral (63) C ηk,ℓ = τk,ℓ

Direct normal-matrix bound Direct normal-matrix bound Uniform symbol bound Normal-matrix bound Resolvent / normal-matrix bound

Exact tanh by quadrature; optional erf surrogate for the synthetic and timing experiments. Same real-valued nonlinear-kernel construction. Same real-valued nonlinear-kernel construction. Same construction with the Fourier coefficients. Same construction with the skew-symmetric coefficients. Linear covariance explicit; nonlinear kernel requires a complex feature convention.

C. Gaussian skew-symmetric matrices Let Gn − G⊤ Sn = √ n , Wr(n) = σr Sn . 2n Since Sn⊤ = −Sn , the Hermitian matrix iSn has the standard semicircular limit [34]. Writing √ 4 − x2 1[−2,2] (x) dx, dµsc (x) = 2π the coefficients are Z 2 skew (1−α+iασr x)k (1−α−iασr x)ℓ dµsc (x). (63) ηk,ℓ = −2

Permutation symmetry and the Gaussian-Poincaré argument of Lemma 2 give diagonal self-averaging. Since An is normal, complete history follows under the asymptotic condition p (1 − α)2 + 4α2 σr2 < 1. D. Complex-valued matrices For complex recurrent matrices, transpose is replaced by the Hermitian adjoint. Normalized complex Ginibre and scaled Haar-unitary matrices again satisfy
1 C Tr Akn (Aℓn )∗ = τk,ℓ ηk,ℓ := lim n→∞ n [34], [37]. The corresponding Hermitian linear-state covariance is obtained from (36) by replacing the input inner product with u∗t−1−k us−1−ℓ . A nonlinear feature kernel additionally requires a specified complex activation or a realification convention. The experiments in this paper use real-valued recurrent matrices. A PPENDIX E F ULL TASK -L EVEL D EPLOYMENT S CORES Table VII gives per-task scores for the four primary selectors, including ESN Direct500 , at all five deployment widths (n ∈ {1 000, 3 000, 5 000, 10 000, 20 000}), while Figure 4 shows the corresponding width trends. A PPENDIX F S ELECTION -L ENGTH S ENSITIVITY AND S ELECTION T IME Table VIII and Figures 5–6 report the auxiliary coarse-grid sweep over the selection-sequence length Tselect summarized in Section IV-B: deployment performance at n = 20 000 and per-task selection CPU time as a function of Tselect , at nselect = 500.

A. Surrogate-kernel validation The main synthetic selection and timing experiments use the closed-form erf surrogate of Remark 1. All real-data forecasting experiments, both ETT and Telco, use the exact tanh kernel evaluated by Gauss–Hermite quadrature. The surrogate is therefore used only as a computational accelerator in the controlled synthetic setting. We validated this choice by comparing matched erfsurrogate and exact-tanh FP selections on the synthetic suite. Across the checked selection lengths, the largest task-level deployment-score difference was 6.3 × 10−3 , and the mean difference was O(10−3 ), well below the task-level variability in the deployment tables. This supports using the surrogate for the synthetic selection and timing sweeps. By contrast, the real-data experiments are reported with the exact tanh kernel, so no real-data conclusion relies on the surrogate approximation. A PPENDIX G F ORECASTING H ORIZON -L EVEL R ESULTS Table IX gives the underlying horizon-level scores. R EFERENCES [1] M. Lukoševičius and H. Jaeger, “Reservoir computing approaches to recurrent neural network training,” Computer Science Review, vol. 3, no. 3, pp. 127–149, 2009. [2] G. Tanaka, T. Yamane, J. B. Héroux, R. Nakane, N. Kanazawa, S. Takeda, H. Numata, D. Nakano, and A. Hirose, “Recent advances in physical reservoir computing: A review,” Neural Networks, vol. 115, pp. 100–123, 2019. [3] H. Jaeger, “The “echo state” approach to analysing and training recurrent neural networks-with an erratum note,” Bonn, Germany: German national research center for information technology gmd technical report, vol. 148, no. 34, p. 13, 2001. [4] H. Zhou, S. Zhang, J. Peng, S. Zhang, J. Li, H. Xiong, and W. Zhang, “Informer: Beyond efficient transformer for long sequence time-series forecasting,” in Proceedings of the AAAI conference on artificial intelligence, vol. 35, no. 12, 2021, pp. 11 106–11 115. [5] J. Bergstra and Y. Bengio, “Random search for hyper-parameter optimization,” Journal of Machine Learning Research, vol. 13, pp. 281–305, 2012. [6] J. Bergstra, R. Bardenet, Y. Bengio, and B. Kégl, “Algorithms for hyperparameter optimization,” in Advances in Neural Information Processing Systems, vol. 24. Curran Associates, Inc., 2011, pp. 2546–2554. [7] J. R. Maat, N. Gianniotis, and P. Protopapas, “Efficient optimization of echo state networks for time series datasets,” in 2018 International Joint Conference on Neural Networks (IJCNN). IEEE, 2018, pp. 1–7. [8] L. A. Thiede and U. Parlitz, “Gradient based hyperparameter optimization in echo state networks,” Neural Networks, vol. 115, pp. 23–29, 2019.

17

TABLE VII F ULL TASK - LEVEL DEPLOYMENT SCORES UNDER THE HOLDOUT- RIDGE PROTOCOL . VALUES ARE 1 − NRMSE MEAN±STD OVER DEPLOYMENT SEEDS . I NTERMEDIATE WIDTHS USE THREE DEPLOYMENT SEEDS ; n = 20 000 USES TEN DEPLOYMENT SEEDS . B EST SCORE PER TASK AND WIDTH IS BOLDED ; SECOND - BEST IS UNDERLINED . I DENTICAL FP AND D IRECT500 ROWS ( E . G ., MC) INDICATE THAT BOTH SELECTORS CHOSE THE SAME CONFIGURATION .

Task MC

Selector

FP Memory500 Direct500 ESN Direct500 NARMA10 FP Memory500 Direct500 ESN Direct500 NARMA20 FP Memory500 Direct500 ESN Direct500 NARMA30 FP Memory500 Direct500 ESN Direct500 NARMA50 FP Memory500 Direct500 ESN Direct500 Inubushi FP Memory500 Direct500 ESN Direct500 Lorenz63 FP Memory500 Direct500 ESN Direct500 SF-MG30 FP Memory500 Direct500 ESN Direct500 MG84 FP Memory500 Direct500 ESN Direct500 Lorenz96 FP Memory500 Direct500 ESN Direct500

n = 1 000

n = 3 000

n = 5 000

n = 10 000

n = 20 000

0.322±0.023 0.418±0.004 0.439±0.001 0.457±0.006 0.463±0.003 0.184±0.004 0.211±0.004 0.214±0.003 0.279±0.005 0.283±0.001 0.322±0.023 0.418±0.004 0.439±0.001 0.457±0.006 0.463±0.003 0.223±0.004 0.268±0.015 0.302±0.003 0.318±0.006 0.328±0.002 0.763±0.014 0.850±0.012 0.854±0.004 0.871±0.010 0.870±0.034 0.720±0.014 0.782±0.008 0.789±0.002 0.804±0.004 0.805±0.030 0.810±0.012 0.820±0.009 0.825±0.009 0.824±0.009 0.832±0.025 0.820±0.011 0.836±0.004 0.865±0.009 0.872±0.007 0.875±0.025 0.738±0.007 0.856±0.010 0.880±0.006 0.909±0.009 0.917±0.009 0.725±0.003 0.824±0.012 0.843±0.006 0.871±0.011 0.877±0.007 0.805±0.005 0.904±0.005 0.915±0.004 0.931±0.007 0.936±0.007 0.739±0.006 0.831±0.008 0.847±0.006 0.862±0.002 0.870±0.009 0.618±0.010 0.730±0.019 0.773±0.008 0.824±0.002 0.835±0.013 0.602±0.022 0.689±0.007 0.734±0.006 0.770±0.006 0.772±0.010 0.613±0.010 0.653±0.056 0.745±0.010 0.790±0.005 0.797±0.011 0.607±0.009 0.608±0.007 0.610±0.008 0.609±0.010 0.605±0.010 0.598±0.005 0.605±0.004 0.606±0.004 0.606±0.005 0.610±0.008 0.412±0.045 0.538±0.006 0.553±0.002 0.574±0.004 0.578±0.012 0.602±0.005 0.606±0.004 0.607±0.004 0.607±0.005 0.610±0.008 0.573±0.008 0.590±0.004 0.596±0.003 0.600±0.004 0.607±0.009 0.948±0.004 0.960±0.003 0.962±0.002 0.962±0.002 0.966±0.003 0.793±0.005 0.798±0.010 0.803±0.007 0.811±0.009 0.819±0.004 0.969±0.002 0.978±0.001 0.979±0.001 0.979±0.002 0.980±0.001 0.994±0.001 0.995±0.001 0.996±0.000 0.995±0.000 0.994±0.000 0.600±0.006 0.620±0.005 0.618±0.003 0.621±0.002 0.624±0.004 0.490±0.013 0.486±0.011 0.496±0.019 0.475±0.003 0.490±0.015 0.652±0.007 0.665±0.004 0.669±0.008 0.672±0.001 0.685±0.007 0.738±0.003 0.733±0.023 0.687±0.012 0.674±0.044 0.679±0.033 0.898±0.014 0.918±0.003 0.931±0.003 0.928±0.002 0.926±0.002 0.841±0.009 0.870±0.020 0.887±0.004 0.891±0.003 0.891±0.004 0.904±0.013 0.925±0.001 0.927±0.004 0.926±0.001 0.927±0.002 0.893±0.001 0.902±0.003 0.884±0.009 0.898±0.003 0.898±0.008 0.812±0.005 0.816±0.001 0.816±0.002 0.819±0.001 0.822±0.002 0.786±0.014 0.802±0.018 0.827±0.003 0.823±0.011 0.828±0.003 0.844±0.018 0.829±0.004 0.833±0.001 0.835±0.000 0.836±0.001 0.827±0.006 0.824±0.005 0.815±0.003 0.816±0.001 0.818±0.007 0.629±0.004 0.666±0.018 0.683±0.001 0.686±0.001 0.686±0.002 0.637±0.017 0.677±0.002 0.680±0.001 0.683±0.001 0.683±0.002 0.679±0.003 0.694±0.001 0.680±0.021 0.666±0.022 0.671±0.023 0.664±0.002 0.668±0.000 0.669±0.000 0.670±0.000 0.670±0.000

[9] F. Matzner, “Hyperparameter tuning in echo state networks,” in Proceedings of the Genetic and Evolutionary Computation Conference, 2022, pp. 404–412. [10] A. Racca and L. Magri, “Robust optimization and validation of echo state networks for learning chaotic dynamics,” Neural Networks, vol. 142, pp. 252–268, 2021. [11] M. Lukoševičius and A. Uselis, “Efficient cross-validation of echo state networks,” in International conference on artificial neural networks. Springer, 2019, pp. 121–133. [12] I. B. Yildiz, H. Jaeger, and S. J. Kiebel, “Re-visiting the echo state property,” Neural networks, vol. 35, pp. 1–9, 2012. [13] M. Lukoševičius, “A practical guide to applying echo state networks,” in Neural Networks: Tricks of the Trade, 2nd ed. Springer, 2012, pp. 659–686. [14] D. Verstraeten and B. Schrauwen, “On the quantification of dynamics in reservoir computing,” in International Conference on Artificial Neural Networks. Springer, 2009, pp. 985–994. [15] B. Schrauwen, M. Wardermann, D. Verstraeten, J. J. Steil, and D. Stroobandt, “Improving reservoirs using intrinsic plasticity,” Neurocomputing, vol. 71, no. 7–9, pp. 1159–1171, 2008. [16] M. Hermans and B. Schrauwen, “Recurrent kernel machines: Computing with infinite echo state networks,” Neural Computation, vol. 24, no. 1, pp. 104–133, 2012. [17] R. Couillet, G. Wainrib, H. T. Ali, and H. Sevi, “A random matrix

approach to echo-state neural networks,” in International Conference on Machine Learning (ICML), 2016, pp. 517–525. [18] J. Dong, R. Ohana, M. Rafayelyan, and F. Krzakala, “Reservoir computing meets recurrent kernels and structured transforms,” in Advances in Neural Information Processing Systems, vol. 33, 2020, pp. 16 785– 16 796. [19] G. A. D’Inverno and J. Dong, “Comparison of reservoir computing topologies using the recurrent kernel approach,” Neurocomputing, vol. 611, p. 128679, 2025. [20] L. Gonon, L. Grigoryeva, and J.-P. Ortega, “Reservoir kernels and volterra series,” IEEE Transactions on Neural Networks and Learning Systems, vol. 37, no. 5, pp. 2181–2192, 2026. [21] A. Gu, K. Goel, and C. Re, “Efficiently Modeling Long Sequences with Structured State Spaces,” in International Conference on Learning Representations, 2022. [Online]. Available: https://openreview.net/ forum?id=uYLFoz1vlAC [22] A. Orvieto, S. L. Smith, A. Gu, A. Fernando, C. Gulcehre, R. Pascanu, and S. De, “Resurrecting recurrent neural networks for long sequences,” in International conference on machine learning. PMLR, 2023, pp. 26 670–26 698. [23] L. Gonon and J.-P. Ortega, “Reservoir computing universality with stochastic inputs,” IEEE Transactions on Neural Networks and Learning Systems, vol. 31, no. 1, pp. 100–112, 2020. [24] L. Grigoryeva and J.-P. Ortega, “Universal discrete-time reservoir

18

Fig. 4. Deployment score vs. reservoir width n through n = 20 000 for the ten-task, four-selector comparison. Shaded bands show ±1 s.d. over deployment seeds; the n = 20 000 column uses ten deployment seeds and the intermediate widths use three. Higher is better. TABLE VIII S ELECTION - LENGTH SENSITIVITY IN AN AUXILIARY COARSE - GRID SYNTHETIC SWEEP. S CORES ARE DEPLOYMENT 1 − NRMSE AT n = 20 000, AVERAGED OVER THE TEN SYNTHETIC TASKS ; ± DENOTES TASK - LEVEL STANDARD DEVIATION . T HE SELECTION LENGTH Tselect IS USED ONLY DURING HYPERPARAMETER SELECTION . FP USES THE CLOSED - FORM ERF - SURROGATE SELECTION KERNEL FROM R EMARK 1; MATCHED EXACT-tanh CHECKS ARE REPORTED IN S EC . F-A. M EMORY 500 USES D = min(Tselect , 1000) DELAYED TARGETS AND SELECTS THE SAME CANDIDATE THROUGHOUT THIS SWEEP. B EST ENTRIES ARE BOLDED AND SECOND - BEST ENTRIES ARE UNDERLINED .

Tselect

FP ↑

Memory500 ↑

Direct500 ↑

ESN Direct500 ↑

150 250 500 1000 2000

0.627±0.196 0.744±0.148 0.771±0.158 0.754±0.185 0.749±0.202

0.727±0.165 0.727±0.165 0.727±0.165 0.727±0.165 0.727±0.165

0.662±0.176 0.688±0.188 0.705±0.147 0.746±0.176 0.761±0.169

0.640±0.261 0.626±0.273 0.698±0.191 0.710±0.190 0.727±0.192

computers with stochastic inputs and linear reservoirs using nonhomogeneous state-affine systems,” Journal of Machine Learning Research, vol. 19, no. 24, pp. 1–40, 2018. [25] G. Lagomarsini, A. Ceni, and C. Gallicchio, “Benchmarking nonlinear readouts in linear reservoir networks,” in International Conference on Artificial Neural Networks. Springer, 2025, pp. 176–187. [26] C. Bordenave and D. Chafaï, “Around the circular law,” Probability Surveys, vol. 9, pp. 1–89, 2012. [27] A. F. Atiya and A. G. Parlos, “New results on recurrent network training: unifying the algorithms and accelerating convergence,” IEEE transactions on neural networks, vol. 11, no. 3, pp. 697–709, 2000. [28] H. Jaeger, “Short term memory in echo state networks,” German National Research Center for Information Technology (GMD), Sankt Augustin, Germany, Tech. Rep. GMD Report 152, 2001. [29] M. Inubushi and K. Yoshimura, “Reservoir computing beyond memorynonlinearity trade-off,” Scientific reports, vol. 7, no. 1, p. 10199, 2017. [30] E. N. Lorenz, “Deterministic non-periodic flow,” J. Atmos. Sci., vol. 20, pp. 130–141, 1963.

[31] M. C. Mackey and L. Glass, “Oscillation and chaos in physiological control systems,” Science, vol. 197, no. 4300, pp. 287–289, 1977. [32] E. N. Lorenz, “Predictability: A problem partly solved,” in Proceedings of the ECMWF Seminar on Predictability, vol. 1, Reading, UK, 1996, pp. 1–18. [33] R. Vershynin, High-Dimensional Probability: An Introduction with Applications in Data Science, ser. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2018, vol. 47. [34] J. A. Mingo and R. Speicher, Free probability and random matrices. Springer, 2017, vol. 35. [35] R. Bhatia, Matrix Analysis, ser. Graduate Texts in Mathematics. Springer, 1997, vol. 169. [36] H. Schultz, “Non-commutative polynomials of independent Gaussian random matrices: The real and symplectic cases,” Probability Theory and Related Fields, vol. 131, no. 2, pp. 261–309, 2005. [37] B. Collins and C. Male, “The strong asymptotic freeness of haar and deterministic matrices,” Annales Scientifiques de l’École Normale Supérieure, vol. 47, no. 1, pp. 147–163, 2014.

19

Fig. 5. Mean deployment score at n = 20 000 as a function of Tselect . Higher is better. FP reaches its best observed mean at Tselect = 500; Direct500 and ESN Direct500 improve with longer sequences. Memory500 selects the same candidate across this sweep and is therefore flat; its proxy uses D = min{Tselect , 1000}.

Fig. 6. Mean selection CPU time per task averaged over the ten-task synthetic suite as a function of Tselect (nselect = 500). All points are measured runs. FP uses the closed-form erf-surrogate selection kernel. Lower is better.

TABLE IX H ORIZON - LEVEL FORECASTING RESULTS AT n = 20 000. S CORES ARE 1 − NRMSE; HIGHER IS BETTER . ETT H IS SAMPLED HOURLY, SO HORIZON STEPS COINCIDE WITH PHYSICAL HOURS ; ETT M IS SAMPLED EVERY 15 MINUTES , AND HORIZONS ARE GIVEN IN STEPS WITH THE PHYSICAL HORIZON IN PARENTHESES . ETT DEVIATIONS ARE OVER THREE DEPLOYMENT SEEDS FOR DETERMINISTIC SELECTORS ; T ELCO DEVIATIONS ARE ACROSS TEN CELLS . T ELCO SELECTION OPTIMIZES THE AGGREGATE OVER HORIZONS 1–12; THE TABLE SHOWS THREE REPRESENTATIVE HORIZONS , SO THEIR MEAN DIFFERS FROM THE AGGREGATE SCORES OF TABLE IV OF THE MAIN PAPER .

Dataset

Metric

FP

Memory500

ESN Direct500

Direct500

ETTh1

1h 6h 12 h

0.695 ± 0.00056 0.597 ± 0.00024 0.731 ± 0.00048 0.724 ± 0.00068 0.587 ± 0.00012 0.576 ± 0.00024 0.586 ± 0.00156 0.589 ± 0.00044 0.386 ± 0.00102 0.392 ± 0.00009 0.374 ± 0.00213 0.387 ± 0.00100

ETTh2

1h 6h 12 h

0.836 ± 0.00022 0.583 ± 0.00010 0.835 ± 0.00063 0.836 ± 0.00022 0.620 ± 0.00135 0.509 ± 0.00005 0.609 ± 0.00165 0.620 ± 0.00135 0.537 ± 0.00027 0.483 ± 0.00013 0.542 ± 0.00144 0.537 ± 0.00027

ETTm1 4 steps (1 h) 0.719 ± 0.00145 0.702 ± 0.00021 0.728 ± 0.00059 0.714 ± 0.00088 24 steps (6 h) 0.446 ± 0.00116 0.435 ± 0.00012 0.449 ± 0.00079 0.438 ± 0.00191 48 steps (12 h) 0.351 ± 0.00123 0.357 ± 0.00003 0.346 ± 0.00097 0.345 ± 0.00089 ETTm2 4 steps (1 h) 0.862 ± 0.00046 0.762 ± 0.00008 0.861 ± 0.00025 0.862 ± 0.00046 24 steps (6 h) 0.456 ± 0.00020 0.474 ± 0.00014 0.457 ± 0.00013 0.456 ± 0.00020 48 steps (12 h) 0.442 ± 0.00004 0.501 ± 0.00025 0.443 ± 0.00012 0.442 ± 0.00004 Telco

1h 6h 12 h

0.601 ± 0.180 0.567 ± 0.162 0.549 ± 0.192

0.679 ± 0.092 0.590 ± 0.130 0.589 ± 0.128

0.644 ± 0.151 0.593 ± 0.134 0.585 ± 0.145

0.685 ± 0.088 0.610 ± 0.119 0.611 ± 0.119

Related documents

• Free-Probability Kernels for Zero-Rollout Hyperparameter Selection in Reservoir Computing#662180
• Free-Probability Kernels for Zero-Rollout Hyperparameter Selection in Reservoir Computing#661130
• Uncertainty Quantification of Next Generation Reservoir Computing with Applications to Memory-Driven Dynamical Systems#1116500
• Structural and dynamical strategies to prevent runaway excitation in reservoir computing.#613484
• Inference of Unknown Dynamical Components Using Next Generation Reservoir Computing: From Chaotic Systems to Climate Data#1051157