Romain Couillet
email: romain.couillet@centralesupelec.

Gilles Wainrib
email: gilles.wainrib@ens.

Tiomoko Hafiz
email: hafiz.tiomokoali@centralesupelec.

Ali

Harry Sevi

A Random Matrix Approach to Echo-State Neural Networks

Recurrent neural networks, especially in their linear version, have provided many qualitative insights on their performance under different configurations. This article provides, through a novel random matrix framework, the quantitative counterpart of these performance results, specifically in the case of echo-state networks. Beyond mere insights, our approach conveys a deeper understanding on the core mechanism under play for both training and testing.

Introduction

Echo-state networks (ESN's) are part of the broader family of recurrent neural networks, specifically dedicated to handling time-series related tasks (such as prediction, non-linear interpolation, etc.) (Jaeger, 2001a;b). Their main feature, as opposed to more conventional neural networks, is to rely on a fixed (but generally randomly chosen) connectivity matrix, the so-called reservoir, and only to enforce network-to-sink edges during the training phase. This reduces overfitting but in turn only allows for shortterm memorization capabilities, unlike backward propagated neural networks that instead target long-term memory.

The fact that the reservoir is chosen once and for all, instead of constantly being updated, eases the theoretical analysis of these networks. As such, by means of a succession of works (e.g., (Jaeger, 2001a;[START_REF] Ganguli | Memory traces in dynamical systems[END_REF]Ganguli Proceedings of the 33 rd International Conference on Machine Learning, New York, NY, USA, 2016. JMLR: W&CP volume 48. Copyright 2016 by the author(s). & Sompolinsky, 2010; [START_REF] Strauss | Design strategies for weight matrices of echo state networks[END_REF][START_REF] Charles | Short-term memory capacity in networks via the restricted isometry property[END_REF]), many qualitative aspects of ESN's have been fairly well understood. In particular, it has been made clear that both the operator norm and the spectral radius of the connectivity matrix play a central role in the ESN performance, that normal versus non-normal connectivity matrices convey strikingly different behavior [START_REF] White | Short-term memory in orthogonal neural networks[END_REF], etc. However, to the best of the authors' knowledge, there has never been an attempt to turn these qualitative considerations into concrete quantitative figures of the end-to-end neural network performance (such as here the reconstruction mean square error).

The objective of the present article is to provide a first theoretical analysis of the (mean-square error) performance of linear ESN's with internal noise for both the training and the testing tasks (an illustration of such a network is depicted in Figure 1). To this end, we shall leverage recent tools from the field of random matrix theory (the applications of which are so far almost not existent in neural networks) and shall conveniently work under the assumption that both the reservoir size n and the training (or test-ing) duration T (or T ) are large and commensurable. The large dimensional framework will induce concentration of measure properties that bring asymptotic determinism in the performance of the random outputs. These results take closed-form expressions which, if not completely explicit, are easily interpreted. Turning the connectivity matrices into (large) random realizations of simple random matrix models, we then further simplify these expressions that turn into elementary formulas.

Among the noteworthy outcomes of these theoretical results, we shall structurally isolate the impact of each system component (normal or non-normal connectivity matrix, independent or not input layer, internal noise level, etc.) on both the ESN training and testing performance. We shall also bridge the qualitative notion of Fisher-memory curve to our quantitative findings in specific scenarios. As an aftermath of our theoretical analysis, we shall finally propose a new practical connectivity matrix design, coined multi-memory matrix that helps handling time series with multiple-scale memory properties.

Note in passing that, beyond their machine learning attractiveness for obvious network stability reasons (although ridge regularized networks without internal noise are often preferred), ESN's with internal noise appropriately model biological short-term memory in the brain. Our results may then also embrace both the neurophysiology and neurocomputation fields.

In the remainder of the article, we shall first introduce the necessary random matrix toolbox used throughout the article (so far not conventional in machine learning), before addressing the performance estimation of network training and testing tasks. Applications of our results will be discussed next, before closing on a reflective discussion. The technical results found throughout the paper are proved in the extended version of the present article [START_REF] Couillet | The asymptotic performance of linear echo state neural networks[END_REF].

The Random Matrix Framework

Before delving into the concrete ESN performance, we shall first consider the elementary objects under study from a random matrix perspective.

We assume here an n-node ESN with connectivity matrix W ∈ R n×n , source-to-reservoir vector m ∈ R n , states x t ∈ R n , t = -∞, . . . , ∞, and internal noise ηε t ∼ N (0, η 2 I n ), fed by a scalar input u t ∈ R. The state evolution equation follows:

x t+1 = W x t + mu t+1 + ηε t+1 .
(1)

The ESN will be trained for a period T and tested for a period T , using a least-square regression approach. Denoting X = [x 0 , . . . , x T -1 ] ∈ R n×T , it shall appear in Section 3 that the mean-square error performance in training relies fundamentally on the matrices

Q γ ≡ 1 T XX T + γI n -1 Qγ ≡ 1 T X T X + γI T -1 (2)
for γ > 0. These matrices, respectively called resolvent and co-resolvent of the Gram matrix 1 T XX T in the operator theory jargon, have been extensively used in random matrix theory for various models of X [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF][START_REF] Pastur | Eigenvalue distribution of large random matrices[END_REF] with multiple applications to engineering notably [START_REF] Couillet | Random Matrix Methods for Wireless Communications[END_REF]. The model of X defined through (1) is not part of those models studied in classical random matrix works, but the technical tools exist to handle it. Precisely, we shall use here the Gaussian framework devised by Pastur [START_REF] Pastur | Eigenvalue distribution of large random matrices[END_REF], based on an integration-by-parts formula for Gaussian random variables and the so-called Nash-Poincaré inequality.

To this end, we first need elementary growth assumptions on the size n and the periods T and T (these, in themselves, are a first result of a "stable" growth regime of ESN's).

Assumption 1 Define the matrix

A = M U with M = [m, W m, . . . , W T -1 m] and U = T -1 2 {u j-i } T -1 i,j=0 . Then, as n → ∞, 1. n/T → c ∈ (0, ∞) and n/ T → ĉ ∈ [0, ∞) 2. lim sup n W < 1 3. lim sup n AA T < ∞ with • the operator norm.
For notational convenience, we define the relation X n ↔ Y n to mean that the (random or deterministic) matrices X n , Y n satisfy a T n (X n -Y n )b n → 0, almost surely, for all deterministic unit norm vectors a n , b n . Under Assumption 1, applying the aforementioned Gaussian framework, we have the following result.

Theorem 1 (Deterministic Equivalent) Let Assumption 1 hold. For γ > 0, and with Q γ and Qγ defined in (2), as n → ∞,

Q γ ↔ Qγ ≡ 1 γ I n + η 2 Rγ + 1 γ A I T + η 2 R γ -1 A T -1 Qγ ↔ Qγ ≡ 1 γ I T + η 2 R γ + 1 γ A T I n + η 2 Rγ -1 A where R γ ∈ R T ×T and Rγ ∈ R n×n are solutions to R γ = 1 T tr S i-j Qγ T i,j=1 Rγ = ∞ q=-∞ 1 T tr J q Qγ S q with [J q ] ij ≡ δ i+q,j , S q ≡ k≥0 W k+(-q) + (W k+q + ) T (and (x) + = max(x, 0)).
Theorem 1 precisely states that any random bilinear form of the type a T n Q γ b n for deterministic a n , b n , can be well approximated, for all large n, T , by the deterministic quantity a T n Qγ b n . As shall be seen in Section 3, this will provide us with a deterministic approximation of the training performance of ESN's. Note in passing that Theorem 1 is strongly reminiscent of the earlier results [START_REF] Hachem | The empirical distribution of the eigenvalues of a Gram matrix with a given variance profile[END_REF] in a close but different context, and are thus not surprising.

The testing phase performance is more involved and does not solely rely on a single Q γ -type matrix. This is because this phase involves both the trained dataset X and the newly observed states X = [x 0 , . . . , x T -1 ], where xt = x t+L for some L T (we assume here a post-training washout step for simplicity). Surprisingly enough, we shall not need the advanced statistics of X explicitly (as opposed to X as seen previously) but shall rather exploit its independence from X. As such, the quantities at stake here are the matrices

1 √ T Q γ X, 1 T X T Q γ BQ γ X
for B any symmetric matrix independent of X such that lim sup n B < ∞. Precisely, we have the following second deterministic equivalents.

Theorem 2 (Second deterministic equivalent) Let Assumption 1 hold and let B ∈ R n×n be a symmetric matrix of bounded spectral norm, independent of X. Then, recalling the notations of Theorem 1, for every γ > 0,

Q γ 1 √ T X ↔ Qγ A(I n + η 2 R γ ) -1 1 T X T Q γ BQ γ X ↔ η 2 γ 2 Qγ G [B] γ Qγ + P T Qγ B + G[B] γ Qγ P with P = A(I n + η 2 R γ ) -1 and where G [B] γ , G[B] γ are so- lutions to G [B] γ = 1 T tr S i-j Qγ B + G[B] γ Qγ T i,j=1 G[B] γ = ∞ q=-∞ η 4 γ 2 1 T tr J q Qγ G [B] γ Qγ S q .
Equipped with these technical results, we are now in position to provide our main contribution to the asymptotic performance of ESN as n, T, T → ∞.

Asymptotic Performance

Recall that the ESN under study is defined by the state equation ( 1). We shall successively discuss the performance of the training and testing steps of this linear ESN.

Training Performance

In the training phase, one wishes to map an input sequence u = [u 0 , . . . , u T -1 ] T to a corresponding known output sequence r = [r 0 , . . . , r T -1 ] T . To this end, we shall enforce the reservoir-to-sink connections of the network, gathered into a vector ω ∈ R n and depicted in color in Figure 1, so to minimize the quadratic reconstruction error

E η (u, r) ≡ 1 T X T ω -r 2 .
The solution to this classical problem is to take ω to be the least-square regressor

ω ≡ (XX T ) -1 Xr , T > n X(X T X) -1 r , T ≤ n. (3)
To such an ω are associated an E η (u; r) which it is convenient to see here as

E η (u, r) = lim γ↓0 γ 1 T r T Qγ r , T > n 0 , T ≤ n (4)
where Qγ was defined in (2).

Applying Theorem 1 in the limit where γ → 0, we have the following first limiting performance result.

Proposition 1 (Training MSE) Let Assumption 1 hold and let r ∈ R T be a vector of Euclidean norm O( √ T ). Then, with E η (u, r) defined in (4), as n → ∞,

E η (u, r) ↔ 1 T r T Qr , c < 1 0 , c > 1.
where, for c < 1,

Q ≡ I T + R + 1 η 2 A T R-1 A -1 and R, R are solutions to 1 R = c 1 n tr S i-j R-1 T i,j=1 R = ∞ q=-∞ 1 T tr J q (I T + R) -1 S q .
Although seemingly not simple, note that, by writing

A T R-1 A = U T m T (W i-1 ) T R-1 W j-1 m T i,j=1
U the matrix Q involved in the asymptotic expression for E η (u, r) clearly features independently:

• the input data matrix U composed in columns of the successive delayed versions of the vector

T -1 2 [u -(T -1) , . . . , u T -1 ] T ;
• the network structuring matrices R and

(W i-1 ) T R-1 W j-1 ;
• the factor η -2 , not present in R, R, which trades off the need for regularizing the ill-conditioned matrix A T R-1 A (through the matrix M in A) and the need to increase the weight of the information-carrying matrix

A T R-1 A (through the matrix U in A).
Note in particular that, since W < 1, the matrix {m T (W i-1 ) T R-1 W j-1 m} T i,j=1 has an exponentially decaying profile down the rows and columns (essentially decaying with i + j). As such, all but the first few columns of R-1 2 M U vanish as n, T grow large, providing us with a first testimony of the ESN short term memory specificity, since only the first columns of U (i.e., the first delays of {u t }) are accounted for. The matrix R-1 2 M then plays the important role of tuning the memory decay.

A particularly interesting scenario is when c = 0 (i.e., n/T → 0). In this case, it is immediate that R = 0 and

R = S 0 = k≥0 W k (W k ) T , so that E η (u, r) ≡ η 2 T r T η 2 I T + U T M T S -1 0 M U -1 r
where

[M T S -1 0 M ] kk = m T (W k-1 ) T S -1 0 W k-1
m is recognized to be J(k -1) with J the Fisher memory curve introduced in [START_REF] Ganguli | Memory traces in dynamical systems[END_REF]. The latter was shown to provide a qualitative measure of the capacity of the ESN to store a k-step delayed information. In Proposition 1, the notion is generalized to account for the finiteness of T with respect to n (i.e., c > 0) and, even when c = 0, conveys a non trivial importance to the off-diagonal terms m T (W i-1 ) T R-1 W j-1 m, i = j.

1 R and R are rigorously the limits of Rγ and γ Rγ from Theorem 1, respectively, as γ ↓ 0.

Testing Performance

Having considered the training performance in Section 3.1, we now assume ω given through its definition (3). To avoid unnecessary complications, we shall stick here to the case where c < 1 and leave the case c > 1 to the extended version of the article [START_REF] Couillet | The asymptotic performance of linear echo state neural networks[END_REF]. We wish to evaluate the performance in testing new data inputs û = [û 0 , . . . , û T -1 ] T which ought to be mapped by the ESN to the desired output r = [r 0 , . . . , r T -1 ] T . Thus our next quantity of interest is the mean square error

Êη (u, r; û, r) = 1 T r -XT ω 2 where X = [x 0 , . . . , x T -1 ] with xt =
x t+L for some L sufficiently large (larger than 2T , say) to ensure approximate independence between {x t } and {x t }. That is, we assume a sufficiently long wash-out period between training and testing (as conventionally done). Alternatively, one may merely reinitialize the network after training and generate a sufficiently long dry-run period prior to testing. Similarly, we shall denote next

 = M Û with Ûij = ûj-i and M = [m, . . . , W T -1 m].
Developing the expression of Êη (u, r; û, r), it is convenient to observe that

Êη (u, r; û, r) = 1 T r 2 + lim γ↓0 1 T 2 T r T X T Q γ X XT Q γ Xr -lim γ↓0 2 T T rT XT Q γ Xr.
We may then straightforwardly apply Theorem 2 in the limit of vanishing γ to retrieve the testing counterpart of Proposition 1 as follows.

Proposition 2 (Testing MSE) Let Assumption 1 hold with c < 1 and let r ∈ R T be a vector of Euclidean norm O( T ). Then, as n → ∞, with the notations of Proposition 1,

Êη (u, r; û, r) ↔ 1 η 2 ÂT QP r √ T - r T 2 + 1 T r T QG Qr + 1 η 2 T r T P T Q S 0 + G QPr with P = A(I T + R) -1 and (G, G) solution to G = c 1 n tr S i-j R-1 S 0 + G R-1 T i,j=1 G = ∞ q=-∞ 1 T tr J q (I T + R) -1 G(I T + R) -1 S q .
Note that the deterministic approximation for Êη (u, r; û, r) in Proposition 2 may be divided into a first term involving û, r and the next two terms only involving u, r. As such, once training is performed, only the former term may alter the asymptotic performance. Again, the case c = 0 leads to simpler expressions, as G = 0 and G = 0 in this case, so that

Êη (u, r; û, r) ↔ ÂT η 2 S 0 + AA T -1 A r √ T - r T 2 + 1 T r T A T η 2 S 0 + AA T -2 Ar.
In order to validate the results of Section 3.1 and Section 3.2, we provide in Figure 2 an example of simulated versus asymptotic performance for a prediction task over the popular Mackey-Glass model [START_REF] Glass | A simple model for phase locking of biological oscillators[END_REF]. Disregarding for the moment the difference between the two displayed theoretical curves, note importantly that the accuracy of the approximation, while increasing as it should for larger values of n, T, T , may strongly decrease as η 2 → 0. This is an expected outcome (from deeper mathematical analysis of the proofs of our propositions) which is reminiscent of the ESN instability observed and documented in (Jaeger, 2001b) as the internal noise vanishes. As a matter of fact, it can be shown that one needs η 2 n -1 2 for a theoretical guarantee that the approximation is accurate. This places our analysis within a regime where overfitting does not occur.

Applications

Let us now move to particular scenarios where Propositions 1-2 either greatly simplify or convey new insights.

Random W matrices

We consider here the scenario where, instead of being considered deterministic, we take W to be a single realization of an elementary (large dimensional) random matrix.

REAL HAAR W .

First consider the scenario where W = σZ with Z random orthogonal with statistical invariance by left-and rightmultiplication by orthogonal matrices, i.e., Z is a random real Haar matrix. Then one can determine an explicit expression for the matrices R, R, G, and G involved in Propositions 1-2. In particular, we have the following result (also depicted in Figure 2).

Corollary 1 (Haar W , c < 1) Let W = σZ with Z random real Haar and m be independent of W with m = 1. Then, under Assumption 1 with c < 1,

E η (u, r) ↔ (1 -c) 1 T r T Qr Êη (u, r; û, r) ↔ 1 η 2 Û T DU Q r √ T - r T 2 + 1 1 -c 1 T r T Qr - 1 T r T Q 2 r where Q = (I T + 1 η 2 U T DU ) -1 , while D ∈ R T ×T and D ∈ R T × T are diagonal with D ii = Dii = (1 -σ 2 )σ 2(i-1) .
We clearly see through Corollary 1 the impact of σ which weighs through D ii the successive delay vectors U •,i starting from i = 1 for zero delay. This is again reminiscent of the works [START_REF] Ganguli | Memory traces in dynamical systems[END_REF] where the diagonal elements of D were understood qualitatively as a memory curve, with the property that i≥1 D ii = 1, so that the ESN allocates a total unit amount of memorization capabilities across the successively delayed versions of u.

This observation inspires the generalization of Corollary 1 to a less obvious, although desirable, structure for W . Indeed, note that, with W = σZ and Z real Haar, memory is allocated according to the exponential decay function k → σ k , thus only allowing for a "single mode" memory, i.e., r t should be an exponentially decaying function of u t-k . If instead, as is more common, r t is a more elaborate function of both close past u t-k but also of further past u t-k values, then it might be appropriate for the ESN not to get restricted to a single k → σ k decay profile.

As such, we propose to consider the matrix W = diag(σ 1 Z 1 , . . . , σ Z ) (with diag the block-diagonal operator), where the matrices Z i ∈ R ni×ni are independent real Haar matrices of given sizes and σ i > 0 assume different values across i. In this case, we have the following natural extension of Corollary 1.

Corollary 2 (Multi-memory W ) Let the assumptions of Corollary 1 hold but for W = diag(σ 1 Z 1 , . . . , σ Z ) with Z i ∈ R ni×ni independent real Haar matrices, i n i = n, and σ i > 0. Then the conclusions of Corollary 1 remain valid but for

D ii = Dii = j=1 n j σ 2(i-1) j j=1 n j (1 -σ 2 j ) -1
.

τ W W + 1 W + 2 W + 3 Figure 3. Memory curve τ → Dτ,τ for W = diag(W1, W2, W3
), Wj = σjZj, Zj ∈ R n j ×n j Haar distributed, σ1 = .99, n1/n = .01, σ2 = .9, n2/n = .1, and σ3 = .5, n3/n = .89. The matrices W + i are defined by

W + i = σiZ + i , with Z + i ∈ R n×n Haar distributed.
10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 -6

= .99Z + 1 , W + 2 = .9Z + 2 , W + 3 = .5Z + 3 (with Z + i
Haar distributed) all defined as in Figure 3, n = 400, T = T = 800. We subsequently move to the next natural model for W , that is W composed of i.i.d. zero mean entries with or without Hermitian symmetry. In the latter case, the study is similar to that of the real Haar case and shall lead to the same result as Corollary 1 but for a different profile of the diagonal entries of the matrix D. As for the former symmetrical case (referred to in random matrix theory as the Wigner case), it leads to a more involved (non explicit) expression for R, which assumes a strikingly different structure than when W is non-symmetric (where R is proportional to the identity matrix). Visually, we find in the large n limit the structure depicted in Figure 5.

When placed in the context of Proposition 1, the observed checkerboard structure for the Wigner case suggests an inappropriate spread of the reservoir energy when it comes to fulfilling pure delay tasks. This is indeed observed numerically with strong performance losses already induced by elementary memorization tasks. A particular example is depicted in Figure 6 where, again for the Mackey-Glass input dataset but for a memorization task (consisting in recalling τ past outputs rather than predicting future outputs). This study and the observed performance results suggest an outstanding performance advantage of non-normal versus normal matrix structures, which might deserve deeper future investigation.

Internal noise and robustness

We complete this section by investigating a scenario in which we assume input data corrupted by extra noise in the testing dataset. This models the fact that one may often possess a large bank of "clean" data to train an ESN but that the reality of test data can sometimes be somewhat, if not strongly, different. Our toy model here consists in considering that an extra Gaussian noise is added to the data û with probability p on each sample ût , while the output r is still expected to be consistent with the noiseless version of û.

From a proper theoretical analysis of the result of Proposition 2, conducted precisely in the extended version of the article, we claim that, letting s 2 be the aforementioned additive noise variance and assuming that Êη → 0 as η → 0 and s = 0 (a scenario that can be enforced by letting e.g., r be a linear combination of finitely many past inputs of u), then there exists a trade-off by which too small values of η 2 induce an increase in the mean square error (all the more that s 2 is large) while large values of η 2 induce too much internal noise. There thus exists an optimal value for η 2 which minimizes the testing MSE.

This phenomenon is depicted in a concrete scenario in Figure 7, still for the same Mackey-Glass one-step ahead prediction task. A particular realization of a random Mackey-Glass time series is also presented in Figure 8, which clearly highlights the robustness strength of internal noise.

Concluding Remarks

The random matrix framework introduced in this article brings new light to the actual performance of linear echostate networks with internal noise (and more generally re- current networks), where existing works so far only provide insights based on partial considerations of the network processing (such as information theoretic metrics of the reservoir information). Our results make it clear what levers should be tuned and optimized upon when designing these networks. Although not presently discussed, an outcome of this study (documented in [START_REF] Couillet | The asymptotic performance of linear echo state neural networks[END_REF]) contradicts some beliefs, such as that suggesting that it is appropriate to take m as one of the leading eigenvectors of W ; we can prove that this choice necessarily leads to poor mean square error performance. But aside from purely theoretical considerations, our results also allow for a fast evaluation of the ESN performance without requiring extensive Monte Carlo simulations which we believe experimenters should find convenient.

We stressed above the words linear and internal noise, as they constitute debatable choices. Regarding internal noise, while being an appropriate model assumption in biological networks, it is often regarded as artificial in machine learning (where a regularized least square ω is chosen to stabilize the network). Since large networks induce concentration of measure phenomena that stabilize the MSE performance of the network, we believe that internal noise (leading to random but equally performing outputs) are sometimes more desirable than deterministic (and thus statically biased) outputs. In the full version of this article, comparisons are performed between both cases. But the utmost limiting aspect of the present work rather lies in the linear character of the state equation (1). It is known, more from experiments and insights than theory, that breaking the linear frontier brings vastly more interesting properties to neural networks, with in particular the possibility for W to exceed one. A necessary next step of our study line is the extension of the random matrix framework to handle the behavior of non linear versions of echo-state networks.

Figure 1 .

 1 Figure 1. Echo-state neural network.

Figure 2 .

 2 Figure 2. Training and testing (normalized) MSE for the Mackey Glass one-step prediction, W defined as in Figure 3, n = 200, T = T = 400 (top) and n = 400, T = T = 800 (bottom). Comparison between Monte Carlo simulations (Monte Carlo) and theory from Propositions 1-2 (Th. (fixed W )) or Corollary 2 (Th. (limit)).

Figure 3

 3 Figure 3 depicts the function k → D kk (which again may be thought of as a memory curve) for the "multimemory" matrix W of Corollary 2 versus elementary random Haar matrices with single σ values. Note the evolution of the function slope which successively embraces the memory curves of each individual Haar matrix composing W . Following the same example, we now compare in Figure 4 the testing performance for the Mackey-Glass onestep prediction task under the multi-memory W versus the composing Haar W +i . It is interesting to see here that the multi-memory matrix is almost uniformly more powerful than the composing matrices and that it matches the performance of the best among the latter. This suggests the possibility to use such a matrix structure in scenarios where

Figure 4 .

 4 Figure 4. Testing (normalized) MSE for the Mackey Glass onestep ahead task, W (multimemory) versus W + 1= .99Z + 1 , W + 2 = .9Z + 2 , W + 3 = .5Z + 3 (with Z + iHaar distributed) all defined as in Figure3, n = 400, T = T = 800.

  4.1.2. NORMAL AND NON NORMAL I.I.D. W .

Figure 5 .

 5 Figure5. Upper 9 × 9 part of R for c = 1/2 and σ = 0.9 for W with i.i.d. zero mean Gaussian entries (left) and W Gaussian Wigner (right). Linear grayscale representation with black being 1 and white being 0.

Figure 6 .

 6 Figure 6. Training (top) and testing (bottom) performance of a τ -delay task for τ ∈ {1, . . . , 4} compared for i.i.d. W versus Wigner W , σ = .9 and n = 200, T = T = 400 in both cases (here on the Mackey-Glass dataset).

Figure 7 .Figure 8 .

 78 Figure 7. Testing (normalized) MSE for the Mackey-Glass onestep ahead task with 1% or 10% impulsive N (0, .01) noise pollution in test data inputs, W Haar with σ = .9, n = 400, T = T = 1000. Circles indicate the NMSE theoretical minima. Error bars indicate one standard deviation of the Monte Carlo simulations.

Acknowledgments

The work of Couillet and Tiomoko Ali is supported by the ANR Project RMT4GRAPH (ANR-14-CE28-0006).


Shankar Mahadevan

Federico Angiolini
email: fangiolini@deis.unibo.it

Michael Storgaard

Rasmus Grøndahl Olsen

Jens Sparsø

Jan Madsen

A Network Traffic Generator Model for Fast Network-on-Chip Simulation

For Systems-on-Chip (SoCs) development, a predominant part of the design time is the simulation time. Performance evaluation and design space exploration of such systems in bit-and cycle-true fashion is becoming prohibitive. We propose a traffic generation (TG) model that provides a fast and effective Network-on-Chip (NoC) development and debugging environment. By capturing the type and the timestamp of communication events at the boundary of an IP core in a reference environment, the TG can subsequently emulate the core's communication behavior in different environments. Access patterns and resource contention in a system are dependent on the interconnect architecture, and our TG is designed to capture the resulting reactiveness. The regenerated traffic, which represents a realistic workload, can thus be used to undertake faster architectural exploration of interconnection alternatives, effectively decoupling simulation of IP cores and of interconnect fabrics. The results with the TG on an AMBA interconnect show a simulation time speedup above a factor of 2 over a complete system simulation, with close to 100% accuracy.

Introduction

An important step in the design of a complex Systemon-Chip is to select the optimal architecture for the on-chip network (NoC). In order to do so, it is imperative to analyze and understand network traffic patterns through simulation. This can be accomplished at various stages in the design flow, from abstract transaction level models (TLM) to bitand cycle-true models. In many cases, only the most detailed models prove capable of capturing important aspects of communication performance, e.g. the latency associated with resource contention. The obvious drawback of these approaches is slower simulation speed.

In this paper, we focus on enabling the exploration of different NoC architectures at the bit-and cycle-true level by We assume as a requirement the availability of a reference SoC design, consisting of IP cores and of a NoC, and of an application partitioned and compiled onto the various IP cores. This application might either be software executing on programmable IP cores, or code synthesized into dedicated hardware. This reference system will be used to collect an initial trace of the IP cores' behavior. In order to increase simulation speed for subsequent design space exploration, we propose to replace the IP cores with Traffic Generators (TG) which emulate their communication at the interface with the network, as illustrated in Figure 1.

The goal is to perform only one reference simulation using bit-and cycle-true simulation models of the IP cores running the target application, and to speed up subsequent variants of that simulation using traffic generators coupled with accurate models of the alternative interconnects only. While the internal processing of IP cores does not need thorough replication by the generators, and can often be modeled by waiting for an amount of cycles between network transactions, the unpredictability of network latency of different NoC architectures may lead to changes in the number and relative ordering of transactions. Thus, traffic generators should have at least some reactive capabilities, as will be explained in Section 3.

In order to capture reactive behavior, we propose a TG implementation as a very simple instruction set processor. Our approach is significantly different from a purely behavioral encapsulation of application code into a simulation device, in analogy with TLM modeling. The TG model we propose is aimed at faithfully replicating traffic patterns generated by a processor running an application, not just by the application; this includes e.g. accurate modeling of cache refills and of latencies between accesses, allowing for cycle-true simulations. At the same time, this approach allows a straightforward path towards deployment of the TG device on a silicon NoC test chip.

To evaluate the TG concept, we have integrated the proposed TG model into MPARM [START_REF] Loghi | Analyzing on-chip communication in a MPSoC environment[END_REF], a homogeneous multiprocessor SoC simulation platform, which provides a bitand cycle-true SoC simulation environment. The current version of MPARM supports several NoC architectures, e.g. AMBA [START_REF] Loghi | Analyzing on-chip communication in a MPSoC environment[END_REF], STBus and the ×pipes [START_REF] Dall'osso | xpipes: A latency insensitive parameterized Network-on-Chip architecture for multi-processor SoCs[END_REF], and leverages ARMv7 processors as IP cores. The use of the OCP [1] protocol at the interfaces between the cores and the interconnect allows for easy exchange of IP cores for TGs, as indicated in Figure 1.

The rest of the paper is organized as follows. Section 2 introduces related work, and is followed by a discussion of the requirements for modelling traffic patterns in Section 3. Section 4 details the TG implementation, and Section 5 describes how communication traces are extracted and turned into programs executing on the TG. Section 6 presents initial simulation results which show the potential of our TG approach. Finally, Section 7 provides conclusions.

Related Work

The use of traffic generators to speed up simulation is not new, and several traffic generator approaches and models have been proposed.

In [START_REF] Lahiri | Evaluation of the traffi c-performance characteristics of System-on-Chip communication architectures[END_REF], a stochastic model is used for NoC exploration. Traffic behavior is statistically represented by means of uniform, Gaussian, or Poisson distributions. Such distributions assume a degree of correlation within the communication transactions which is unlikely in a SoC environment. Traffic patterns in SoC systems have shown to be reactive and bursty [START_REF] Bolotin | QNoC: QoS architecture and design process for network on chip[END_REF][START_REF] Lahiri | Communication architecture tuners: A methodology for the design of high-performance communication architectures for System-on-Chips[END_REF]. The simplicity and simulation speed of stochastic models may make them valuable during preliminary stages of NoC development, but, since the characteristics (functionality and timing) of the IP core are not captured, such models are unreliable for optimizing NoC features.

A modeling technique which adds functional accuracy and causality is transaction-level modeling (TLM), which has been widely used for NoC and SoC design [START_REF] Fummi | A timing-accurate modeling and simulation environment for networked embedded systems[END_REF][START_REF] Grötker | System Design with SystemC[END_REF][START_REF] Ogawa | A practical approach for bus architecture optimization at transaction level[END_REF][START_REF] Pasricha | Extending the transaction level modeling approach for fast communication architecture exploration[END_REF][START_REF] Sgroi | Addressing the System-on-Chip interconnect woes through communicationbased design[END_REF]. In [START_REF] Ogawa | A practical approach for bus architecture optimization at transaction level[END_REF][START_REF] Pasricha | Extending the transaction level modeling approach for fast communication architecture exploration[END_REF], TLM has been used for bus architecture exploration. The communication is modeled as read and write transactions which are implemented within the bus model. Depending on the required accuracy of the simulation re-sults, timing information such as bus arbitration delay is annotated within the bus model. In [START_REF] Pasricha | Extending the transaction level modeling approach for fast communication architecture exploration[END_REF] an additional layer called "cycle count accurate at transaction boundary" is presented. Here, the transactions are issued at the same cycle as that observed in bus-cycle-accurate models, thus intratransaction visibility is traded-off for simulation speedup. While modeling the entire system at higher abstraction i.e. TLM, both [START_REF] Ogawa | A practical approach for bus architecture optimization at transaction level[END_REF] and [START_REF] Pasricha | Extending the transaction level modeling approach for fast communication architecture exploration[END_REF] present a methodology for preserving accuracy with gain in simulation speed.

We would like to underline that our approach is dual with respect to TLM. While transaction-level models usually represent interconnects as a collection of available services and emphasize local processing on IP cores, the platform we describe is composed of accurate models for the interconnect, while processing resources are abstracted away. Simulation speed is gained like in TLM models, but the purpose of this gain is enabling accurate assessment of interconnect performance, not of core or application performance. The above methods are suitable for feature exploration once the NoC architecture has been chosen, but are not thought for NoC exploration itself.

Traffic Modeling Requirements

The generation of a traffic pattern emulating that of a real IP core can be faced at varying degrees of accuracy.

At the most basic level, a trace with timestamps can be collected in the reference system and then be independently replayed, an approach that we might call "cloning". This approach is clearly inadequate when the variance of network latency is taken into account; whenever a transaction is delayed, either due to hardware design or congestion, the effect should propagate to subsequent transactions, which would also be delayed in real systems. A simple example of such critical blocking is a cache refill request.

This observation leads to the deployment of "timeshifting" traffic generators: adjacent transactions are tied to each other, and are issued at times which are a function of the delay elapsed before receiving responses to previous transactions. This implicitly means that the trace collection mechanism must include not only timestamps for processor-generated commands, but also for network responses. However, even this model fails when multi-core systems come under scrutiny: the arbitration for resources in such designs is timing-, and thus architecture-, dependent. Therefore, very different transaction patterns may be observed as a function of the chosen interconnection design. To make an example, checks for a shared resource done by polling generate different amounts of traffic depending on the relative ordering of accesses to the resource.

As a consequence, the need for "reactive" TG models is justified. Such models must have some knowledge about the system architecture and about the application behavior to correctly generate (and not just duplicate) traffic patterns across different underlying networks. A TG should be able to mimic the behavior of an IP core even when facing unpredictable network performance, e.g. due to resource contention, packet collisions, arbitration and routing policies.

To illustrate the requirements driving the development of our TG model, we will now describe how two typical transactions occurring in the MPARM modelling and simulation environment can be reproduced. MPARM features in-order, single-pipeline ARM cores as system masters and two types of memory as system slaves: private (only accessible by one master) or shared (visible by all masters in the system). Figure 2 shows examples of the two types of communication: (a) processor-initiated communication towards an exclusively owned slave peripheral, and (b) processorinitiated communication towards a system-shared slave peripheral. We call network latency (t nwk ) the time taken for the communication to traverse from the master OCP interface to the slave OCP interface and vice versa; this latency depends both on the chosen architecture and on the network congestion at the time of communication.

In Figure 2(a), the first two master transactions are a write (WR) and a read (RD). The time to service the WR transaction, which is a posted write, is just the network latency plus the slave access time. The RD, which in MPARM uses blocking semantics, pays an additional penalty because the response has to make its way back to the master. From the TG point of view, this pattern is easily recordable: network latency and slave access time are unimportant factors, and the essential point to capture is just the delay between WR assertion and RD assertion, and between RD response and the following command. In a subsequent simulation with traffic generators replacing cores, these delays will be modeled by explicit idle waits in the TG, while the network latency will be dependent on the NoC model to simulate. In the next set of transactions, where a RD closely follows a WR, the RD command reaches the slave before the latter has finished servicing the WR, and is thus stalled at the slave interface. This stalling behavior does not need to be explicitly captured in a TG model, since, from a processor perspective, it simply appears to be part of the slave response time. This simplistic example of a master accessing a private slave proves that if the type and the timestamp of the communication events are captured, the behaviour of the master can be emulated via non-preemptive sequential communication transactions interleaved with an appropriate amount of idle wait cycles.

In Figure 2(b), two master devices (M1 and M2) attempt to gain access to a single hardware semaphore. M1 arrives first and locks the resource; the attempt by M2 thus fails. In MPARM, semaphore checking is performed by polling, i.e. M2 regularly issues read events until eventually the semaphore is granted to it. Since the transactions occur over a shared network fabric, the unlock event (WR) issued by M1 and the success of the next request (RD) event by M2 are dependent. Only if the M2 RD event is issued at least t nwk,M1 + t unlock,S -t nwk,M2 after the unlocking by M1, then M2 will be granted the semaphore and additional polling events will not be required. Therefore, depending on network properties, a variable amount of transactions might be observed at the OCP interfaces of M1 and M2. This is the reactive behavior that needs to be captured by the TG model: both M1 and M2 need to react to accommodate the network latency. Thus, the simplistic model which could be applied to transactions towards a privately owned slave now needs to be extended with additional information about the master process execution, about system properties, and about input/output data. In detail, the TG must be able to recognize polling accesses (i.e. a knowledge of what addressing ranges represent pollable resources) and must add support for recording of actual data transfers (e.g., writing a "1" or a "0" to a shared memory location might be the difference between locking or unlocking a resource).

We take the above discussion as a requirement to implement accurate TG models. The examples in Figure 2 demonstrate that traces collected at the IP-NoC interface are sufficient to accurately reproduce the IP's communication, provided that the reactive behavior of the master IP cores is taken into account. These traces should collect sequences of communication transactions, comprising of requests and responses, separated by time intervals with no communication, i.e. idle time. A simulation of the entire system should produce several traces, one per IP core interface.

Implementation of the traffic generators

In this section we describe in some detail the implementation of our traffic generators. As mentioned before, our proposed TG model is designed as a simulation tool, but allows future deployment as a hardware device. Within the simulation environment for NoC exploration, the emphasis is on simulation speedup, while within a hardware instance the emphasis is on ease of (re)programmability and a small silicon footprint in order to support implementation of test chips containing NOC prototypes.

Conceptually, three different TG entities might be needed: (1) A TG emulating a processor (an OCP master). This TG must be able to issue conditional sequences of traces composed of communication transactions separated by idle/wait-periods; (2) A TG emulating a shared memory (an OCP slave). This TG must contain a data structure modeling an actual shared memory (since the values read by the masters may affect the sequence of transactions seen at the master IP cores); and (3) A TG emulating a slave memory (an OCP slave). This TG must be able to respond, possibly with dummy values, to communication transactions issued by a master. Only the first is actually required for deployment in a simulation environment, which already provides its own system slaves, thus only this entity will be described in the present paper. On the other hand, it is important to notice that both slave TG modules are much simpler in design with respect to the master TG, as their logic basically just involves a small state machine to handle OCP transactions.

For the processor TG, we have implemented and modeled a multi-cycle processor with a very simple instruction set as listed in Table 1. The processor has an instruction memory and a register file, but no data memory. The instruction set consists of a group of instructions which issue OCP transactions (whose arguments are set up in registers) and a group of instructions allowing the programming of conditional sequencing and parameterized waits such that the required traces can be implemented/programmed. The process for deriving TG programs from traces obtained in the reference simulation is explained in Section 5.

The TG simulation flow

In order to use the traffic generators, a user must first perform a reference simulation using bit-true and cycle-true IP models. It is interesting to note that, at this stage, the interconnect does not yet need to be accurately modeled, allowing for time savings. During this simulation, traces are collected from all OCP interfaces in the system. For this purpose, the OCP interface modules within the MPARM platform (the network interfaces in the case of the ×pipes interconnect, the bus master in the case of AMBA AHB) were adapted to collect traces of OCP request and response communication events into a predefined file format (.trc). The address and (if any) data fields of the transactions are also observed. Trace entries are single or burst read/write transactions. Figure 3 loaded into the TG instruction memory and executed. Execution might be within a simulation model (which is the approach presented in this paper) or in hardware on a NoC test-chip. Validation of the trace collection and processing mechanism can be achieved by collecting traces with IP cores running on different interconnects, and verifying the resulting .tgp and .bin programs to match. The conversion process is fully automated and the time taken for this process is discussed in Section 6. As seen in Figure 3(b), the TG program starts with a header describing the type of core and its identifier. The next few statements express initialization of the register file. Register rdreg is defined as special register where the value of RD transactions is stored.

By looking at the code in Figure 3(a), it is possible to notice that the first communication events in the trace occur at time 55ns, 75ns, and 90ns. We assume each TG cycle to take 5ns, the same as the IP core for which the trace is collected. At the beginning of the simulation, the TG has no instruction to perform until the 11th (55/5) cycle, so an Idle wait is observed. The trace of the RD event is followed by a response, at a time which is dependent on the network latency. The IP core is blocked until this response arrives. A WR event occurs three ((90-75)/5) cycles after the response is received; these cycles are partially spent for TG internal operations (data and address register setting), and an ensuing Idle wait is added to fill the gap. Then the following RD instruction is translated into the corresponding Read program call after 10 cycles, one of which is taken to set up the RD address. This is blocking until a response is received, 5 cycles later. Now, consider the trace entries from time 210ns to 320ns. By identifying the address as belonging to a semaphore location and knowing the polling behaviour of the MPARM IP core, the translator inserts the Semchk label and an If conditional statement. This statement checks whether the read value is equal to "1", which reflects an unblocked semaphore. This loop effectively represents the semaphore polling behavior. All master devices attempting to access this address incorporate the same routine in their TG program, thus capturing the system dynamics.

At this stage, additional simulations can be run on a plat-; Simple RD/WR/WRNP RD 0x00000104 @55ns Resp Data 0x088000f0 @75ns WR 0x00000020 0x00000111 @90ns RD 0x00000031 @140ns Resp Data 0x00002236 @165ns .. .. ; polling a semaphore!! RD 0x000000ff @210ns Resp Data 0x00000000 @270ns RD 0x000000ff @285ns Resp Data 0x00000000 @310ns RD 0x000000ff @305ns Resp Data 0x00000001 @320ns ..

Network latency

Next IP comm transaction interval (a) form with traffic generators and a variety of interconnect fabrics, thereby evaluating performance of NoC design alternatives. Compared with the reference setup, where the interconnect fabric could be modeled at a high level, the target NoC should now be simulated at the cycle-and bittrue level to carefully assess its performance. Validation of the TG model can be achieved by coupling the TG with the same interconnect used for tracing with IP cores, and checking the accuracy of the IP core emulation. Results for this validation, and for tests on different interconnects than the reference one, will be presented in the next Section.

Results

We simulated within the MPARM framework, using the AMBA NoC, and four benchmarks. The first benchmark was a single-processor application (SP matrix manipulation), with the purpose of assessing accuracy and speedup in the simplest environment. The second benchmark (Cacheloop) was a test performing idle loops within the processors' cache, and only minimal bus interaction; this allowed an assessment of the speedup provided by the TG model when scaling the number of processors in the system up to twelve. Finally, the remaining two benchmarks (MP matrix manipulation and DES encryption/decryption) were multiprocessor tests stressing synchronization and resource contention with traffic patterns as discussed in Section 3, and were used mainly to ascertain the accuracy of the whole design flow when stressed by complex transactions.

In the first experiment we aimed at validating the trace collection/processing environment. We ran the same benchmarks over AMBA and ×pipes, noticing very different ex- ecution times due to different latency and scalability features. However, after translation, a check across .tgp programs showed no difference at all. This result demonstrates the feasibility of an approach which decouples simulation of the IP cores and of the underlying interconnect fabric. Table 2 summarizes the results of simulations done on the AMBA AHB interconnect with ARM processors and then with TGs. The left columns report the number of simulated cycles, while the right ones illustrate simulation time 1 . The column "Error" is a measure of the accuracy of replacing IP cores with TGs, based upon the difference in simulated cycles, while the column "Gain" describes the improvement in simulation time.

The table shows that replacing ARM processors with TGs yields excellent accuracy, close to 100% for small numbers of processors, while guaranteeing a speedup factor of 2 to 4. This speedup is mostly due to the drastic simplification in the amount of logic needed to generate communication transactions, compounded in small part with the elimination of any adaptation layer in the system since the TG is natively implemented with an OCP interface. This speedup compares favorably to previous work in the area (a speedup of 1.55x is reported in [START_REF] Pasricha | Extending the transaction level modeling approach for fast communication architecture exploration[END_REF]), and must be evaluated by taking into account the fact that it involves no shift in the level of abstraction of the simulations.

Inaccuracies in execution time can be explained as follows. In Cacheloop, and until about 6-8 processors in MP matrix and four in DES, the TG platform shows an acceptable accuracy degradation. This is due to the compounding of minimal timing mismatches caused by the conversion from traces to TG programs. When adding even more processors, however, accuracy improves again, because the AMBA bus starts to saturate, causing the processors to idle wait for bus arbitration for long amounts of time. The congestion is serious enough to dominate the effect of the timing mismatches. Since TGs cannot save simulation complexity if the replaced processors are in idle state, this is also the reason causing the speedup to get smaller with large numbers of processors in MP matrix and DES. For Cacheloop, which always executes from the local caches without any bus traffic, this phenomenon does not appear. Thus, the reduced speedup is not a property of the TG.

The impact of trace collection is small, and is incurred only once. For example, when running the MP matrix benchmark on the AMBA interconnect with four ARM processors, a plain benchmark run takes 128 s; the benchmark run with TG tracing enabled takes 147 s, and subsequent parsing and elaboration requires an additional 145 s for a 20 MB trace file 1 . Only one such iteration is needed to be able to take advantage of 2x to 4x speedups in subsequent design space exploration. Additionally, since processed TG programs are identical regardless of the reference interconnect in which raw traces were collected, such collection could be performed on top of a transactional fabric model, further reducing the impact of the reference simulation.

Conclusions

Experimental results prove the viability of a TG-based approach which decouples simulation of IP cores and of interconnect fabrics. Even in presence of unpredictable contention for shared resources in a multiprocessor environment, our TG model proved capable of delivering speedups in the order of 2x to 4x when run on AMBA while keeping a remarkable accuracy.

The TG model we propose provides a wide range of features, with a simple but powerful instruction set allowing for sophisticated flow control and therefore a variety of communication patterns. It is very useful for fast and accurate verification and exploration of different NoC architectures, which is the motivation of this work. While this paper was focused on simulation speedup, the TG may also be used as a flexible tool in a variety of platforms. The TG might be used in association with manually written programs to generate traffic patterns typical of IP cores still in the design phase, helping in the tuning of the communication performance between the underlying NoC and that IP core.

Our future work includes synthesis of the TG device, and support for processors allowing out-of-order transac-tions. Research will also include analysis of the behavior of a system in which multiple tasks run on a single processor and are dynamically scheduled by an OS, either based upon timeslices (preemptive multitasking) or upon transition to a sleep state followed by awakening on interrupt receipt. Context switching-related issues will need to be modeled or predicted.

Figure 1 .

 1 Figure 1. Simulation Environment with bitand cycle-true: (a) IP-cores, (b) TG model. increasing the speed of the complete SoC simulation. This is a key advantage, since architectural exploration typically involves carrying out the same set of simulations for each design alternative, and simulations may consist of millions of clock cycles each.We assume as a requirement the availability of a reference SoC design, consisting of IP cores and of a NoC, and of an application partitioned and compiled onto the various IP cores. This application might either be software executing on programmable IP cores, or code synthesized into dedicated hardware. This reference system will be used to collect an initial trace of the IP cores' behavior. In order to increase simulation speed for subsequent design space exploration, we propose to replace the IP cores with Traffic Generators (TG) which emulate their communication at the interface with the network, as illustrated in Figure1.The goal is to perform only one reference simulation using bit-and cycle-true simulation models of the IP cores running the target application, and to speed up subsequent variants of that simulation using traffic generators coupled with accurate models of the alternative interconnects only. While the internal processing of IP cores does not need thorough replication by the generators, and can often be modeled by waiting for an amount of cycles between network transactions, the unpredictability of network latency of different NoC architectures may lead to changes in the number and relative ordering of transactions. Thus, traffic generators should have at least some reactive capabilities, as will be explained in Section 3.

Figure 2 .

 2 Figure 2. Two typical MPARM transactions.

  (a) shows an example trace. The next step is to convert the traces into corresponding TG programs (.tgp). A translator outputs symbolic code; Figure 3(b) shows the TG program derived for traces in Figure 3(a). Finally, an assembler is used to convert the symbolic TG program into a binary image (.bin) which can be

Figure 3 .

 3 Figure 3. (a) MPARM Trace to (b) TG Program.

Table 1 . OCP-master TG instruction set.

 1

	Instructions	Description
	OCP Instructions:
	Read(addr)	Read from an address
	Write(addr, data)	Write to an address
	BurstRead(addr, count)	Burst read a range of addr.
	BurstWrite(addr, data, count)	Burst write an address set
	Other Instructions:
	If(arg1, arg2, operand)	Branch on condition
	Jump(location)	Branch direct
	SetRegister(reg, value)	Set register (load immediate)
	Idle(counter)	Wait for given no of cycles

Table 2 . TG vs. ARM performance with AMBA.

 2

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

Benchmarks taken on a multiprocessor Xeon ® 1.5 GHz with 12 GB of RAM, eliminating any disk swapping effect. Especially for benchmarks with a short duration, time measurements were taken by averaging over multiple runs and care was put in minimizing disk loading effects.Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

Acknowledgments

The work of Shankar Mahadevan is partially funded by SoC-Mobinet, Nokia, and the Thomas B.Thrige Foundation. The work of Federico Angiolini is partially funded by ARTIST.


Stéphanie Fryars

Emmanuelle Limanton

Fabienne Gauffre

Ludovic Paquin

Corinne Lagrost

Philippe Hapiot

Rodger Parsons

DIFFUSION OF REDOX ACTIVE MOLECULES IN DEEP EUTECTIC SOLVENTS

Keywords: Deep Eutectic Solvents, Redox Probes, Cyclic Voltammetry, Diffusion Coefficients

Classical electrochemical approaches were evaluated for investigating the basic properties of deep eutectic solvents (DES) with a special focus on the mass transport phenomena and diffusion of charged species in DES. Three classical DES were considered: Reline, Ethaline and cholinechloride-lactic acid mixture. Cyclic voltammetries of the ferrocene oxidation and of the ferricyanide reduction, respectively with charge 0 and -3, exhibit well-defined responses in these media. For all considered DES and redox molecules, diffusion processes govern the masstransport of the redox molecules from the solution to the electrode. Derived diffusion coefficients were extracted and found to present some inconstancies with the classical Stokes-Einstein law for homogeneous media. The diffusion coefficients were found to be almost unaffected by the charge present on the diffusing molecules contrarily to observations made for ionic liquids.

Introduction

Deep eutectic solvents (DESs) represent a novel class of solvents that has already been proposed for many applications in chemistry including uses in electrochemistry. [1,2,3,4,5] Historically, researches about DESs are related to the Ionic liquids (ILs) that are organic salts with a low melting point (<100°C) combining an organic cation with an organic or inorganic anion. [6,7,8,9,10] DESs and ILs have both attractive physical properties: they are non-volatile, non-flammable and their solvent properties could be tailored by designed functions, therefore the considered applications of DESs or ionic liquids are often similar. They are well-adapted for numerous applications ranging from catalysis, biocatalysis, electrochemistry for energy storage to extraction processes. [6][7][8][9][10][11]12] In a pioneering work, A.P. Abbott et al. have reported the first associations between choline chloride and natural compounds such as urea in 2003. [3] A C C E P T E D M A N U S C R I P T

Associations in DESs could involve different interactions such as hydrogen bonds between a hydrogen-bonding donor (HBD) and a hydrogen-bonding acceptor (HBA). The choice of the components to form a DES is thus immense as DES could associate one to several compounds in different proportions. [1][2][3][4][5]13,14,15,16,17]]

Electrochemistry studies and more particularly cyclic voltammetry have largely been used to characterize the basic properties of neoteric solvents as ionic liquids [7,8] approaches based on cyclic voltammetry and SECM. [13] More recently, the electrochemistry of ferrocene and cobaltocene were investigated using cyclic voltammetry in ChCl based DESs and their diffusion coefficients were derived in these media. [18,19,20] The electrochemical behavior of metal couples such as Fe 3+ /Fe 2+ [ [21,22,23] or Cu 2+ /Cu + [24] were also examined in DES.

They were found to form complexes in ChCl-based DESs making their use as "inert" redox probe more complicated because the diffusive species could be totally different when changing the DES.

Based on NMR investigations, self-diffusion was examined in four ChCl-based DESs and diffusion coefficients were measured. [25] From comparison of experiments with viscosity data and models, it was concluded that the diffusion mechanism is similar to that in pure ionic liquids, which acts via a hopping mechanism or hole diffusion rather than the random walk model used to describe diffusion in conventional homogeneous liquids. [25] In this context, we have performed a series of experiments with a special focus on the diffusion processes involved in DESs. A central question concerns the effect of the charge carried by the redox probe and the inner sphere character of the electron transfer of the redox probe on the mass transport. One could expect that a charged molecule will present a specific association with the cation/anion of the DES providing information about its structure if any. [7,8] We have used the cyclic voltammetry on a millimetric electrode and have selected the reduction of ferricyanide by comparison with the oxidation of ferrocene in 3 different DESs prepared from mixture of choline chloride (ChCl) and different hydrogen bond donors (HBDs) (See Table 1). These two redox molecules were chosen for their chemical stability, the specificity of their electron transfers and the charges present on the molecules. The ferrocene/ferrocenium couple is an outsphere system that is relatively insensitive to solvation

A C C E P T E D M A N U S C R I P T

3 change and its oxidation involves the passage from a charge 0 to +1. In contrast, the reduction of ferricyanide [Fe(CN) 6 ] 3-to ferrocyanide [Fe(CN) 6 ] 4-involves an electron transfer with a strong inner sphere character for a molecule that is negatively charged (passage from -3 to -4 upon reduction). We present here our first results using cyclic voltammetry and preliminary conclusions about the specificity of the mass transport/diffusion properties in DES notably for neutral and charged redox probes.

Experimental.

All electrochemical measurements were performed with an Autolab PGSTAT 302N (Metrohm) using a conventional 3-electrode setup. [7,8] The reference electrode was a quasireference electrode made with a Pt wire. The counter-electrode was a large Pt wire and the working electrode a 1-millimeter diameter Pt disk electrode. Potentials are reported versus the Pt quasi-reference electrode scale for an easier comparison of original experiments.

Maximum faradaic peak current, I p were measured after correction of the baseline that was estimated under the faradaic peak by a linear extrapolation of the background current determined at the beginning of the voltammogram. D values were derived assuming a reversible electron transfer. [START_REF] Savéant | Elements of Molecular and Biomolecular Electrochemistry[END_REF] Ferrocene and potassium ferricyanide were of the highest available purity grade from commercial source (Aldrich) and used without further purification. Solutions of ferrocene in the DESs were prepared by heating the DES with the solute with a gradual increase of the temperature up to 60°C in an ultrasonic bath. We found that concentrations around 5 x 10 -3 mol L -1 are accessible without the appearance of a precipitate when the solution was cooled down to room temperature. Such concentrations were sufficient to record currents with a good ratio between the signal and background currents and noise.

A C C E P T E D

M A N U S C R I P T 4
Ethaline, Reline and choline chloridelactic acid were prepared from commercially available compounds obtained from Aldrich according to published procedures. [START_REF] Hayyan | [END_REF]28] Mass fraction of residual water was measured by using a coulometric Karl-Fisher titration (831KF Coulometer -Metrohm and using Hydranal® Coulomat E solution from Fluka) for each sample and found in the order of 3-7 wt % depending on the DES.

Dynamic viscosities of DES were determined for each sample and the considered temperatures using an Anton-Parr Rheometer MCR301. Experiments are performed with a plateplate setup (diameter: 75 mm and gap: 1mm) in a steady-shear mode. In this protocol, the shear stress and the viscosity are measured as a function of the imposed shear rate. The temperatures is controlled by a Peltier device, with a 0.1°C accuracy.

Results and Discussion

In this work, we use a simple electrochemical system, based on millimetric electrodes and without using electronic compensation of the ohmic drop. Cyclic voltammetries were recorded at different scan rates and at increasing temperature to get a characterization of the mass transport processes in these DESs. DESs were used as prepared to be as close as possible to the solvents used in practical applications. They contain considerable amount of water (3 to 7 wt% depending on the DES). As reported for ILs, the properties of "wet" DES has been shown to vary in a large extent from ultra-dry liquids.

[29]

Investigations in Ethaline (ChCl/ethylene glycol 1:2).

Cyclic voltammetries of the reduction of ferricyanide and oxidation of ferrocene in Ethaline are displayed in Figure 1 at increasing scan rates from 0.02 to 0.5 V.s -1 and for two temperatures (20°C and 60°C). No additional supporting electrolytes were added in the solution meaning that the conductivity of the DES relies solely on its intrinsic conducting properties. It is remarkable that both ferrocene oxidation and ferricyanide reduction display well-defined voltamograms with the classical shape expected for a reversible redox couple. [START_REF] Savéant | Elements of Molecular and Biomolecular Electrochemistry[END_REF] Variations of the faradaic peak currents intensities corrected from the baseline, I p , were measured as function of the scan rate v, (See Figure 2). We used the definition of the peak current Ip for a reversible voltammogram (I p = 0.446 FSC° (Fv/RT) 1/2 to derive the value of the diffusion coefficient; where S is the electrode area, C° the initial concentration of ferrocene or ferricyanide, T the temperature and F and R the Faraday and molar gas constants). [START_REF] Savéant | Elements of Molecular and Biomolecular Electrochemistry[END_REF] For both ferricyanide reduction and ferrocene oxidation, I p varies linearly with v 1/2 between 0.02 -1 V.s -1 range and for the experiments performed at 20°C and 60°C (see Figure 2). Such linear variations mean that the mass-transport of the redox species from the solution bulk to the electrode surface is controlled by their diffusion. From the linear slopes of the I p /v 1/2 variations, we derived the values of the diffusion coefficients, D, of ferricyanide and ferrocene in the DES (See Table 2). D values around 10 -7 cm 2 s -1 are obtained at 20°C and in the 3-4 x 10 -7 range at 60°C which correspond to an increase of the diffusion coefficient by a factor 3.5-3.7 between 20 and 60°C. It is remarkable that the D values for ferricyanide and ferrocene are similar despite the large differences of charges of the two molecules (-3 for ferricyanide and 0 for ferrocene) and display similar increases with the temperature.

Viscosity and amount of water present in the DES (measured on the same sample) are also reported in Table 2. By comparison with literature, derived diffusion coefficient of ferrocene in Ethaline are 2-3 times larger than those previously reported in the same DES at 25°C. [18] These difference falls in line with a lower measured viscosity of our DES; 44 mPas s versus 66 mPas s [18] that could be ascribed to a higher quantity of water in the DES. As done previously (see for example ref 18, 29, 30 and references therein), comparison between D values and the

A C C E P T E D M A N U S C R I P T

6 experimental viscosities brings valuable information about the diffusion mechanism. For an ideal liquid, the classical Stokes-Einstein relation predicts that D varies as T/(where  is the dynamic viscosity of the DES). The validity of the law could be evaluated by simply considering the correlation coefficient =T/( D) that must remain constant with the temperature. [28,30] In ultra-dry ILs and DES, inconstancies versus the Stokes-Einstein law were reported for the diffusion of ferrocene. They were explained by considering the inhomogeneity of such media versus a classical solvent. [25,30] We did not evidence such a discrepancy with our measurements in Ethaline (See Table 2),  varies by around 25 % that is a reasonable agreement with the Stokes-Einstein law considering the uncertainties on the viscosity measurements and diffusion coefficients. At this point, it is also relevant to compare our "electrochemical" data with the self-diffusion coefficients measured by NMR in the same DES and temperature. [25] D value of 4.8 x 10 -7 cm 2 s -1 was reported for the self-diffusion of choline cation at 60°C that is comparable to the 2.3 x 10 - 7 cm 2 s -1 and 3.1 x 10 -7 cm 2 s -1 that we derived for ferricyanide and ferrocene respectively at the same temperature in Ethaline. It is also noticeable that the electrochemical and NMR analyses are performed with different characteristic times, they are typically in the 0.2 -2 s range in our electrochemistry experiments that is longer than the 0.05 s delay used in NMR measurements.

[25]

Investigations in Reline (ChCl/urea 1:2).

Similar electrochemical studies of the reduction of ferricyanide and the oxidation of ferrocene were performed in Reline (see Figure 3) at increasing scan rates from 0.02 to 0.5 V.s -1 . Reline is a less adapted solvent as it presents a much higher viscosity (835 mPa s) at 25 °C and 231 mPa s at 60°C with 3 wt% of water) and a lower conductivity. Additionally, ferrocene was found to As seen in Figures 3A and3B, well-defined voltamograms are obtained for ferricyanide reduction but less-defined curves with large background current for ferrocene oxidation. We could ascribe this observation to partial blockings of the electrode related to the low solubility of ferrocene that favors the deposition onto the electrode surface. At this point, only the electrochemical responses of ferricyanide reduction were analyzed in Reline.

Figure 4 shows the variation of the peak currents with the square root of the scan rate (v 1/2 ).

As for the experiments performed in Ethaline, the peak current displays a clear linear variation with the square root of the scan rate at both 20°C and 60°C. From the slopes, we derived D values of 6.2 x 10 -9 and 5.4 x 10 -9 cm 2 s -1 for ferricyanide in Reline. These values are much lower than those measured in Ethaline as expected from the higher viscosity of Reline. However, the correlation coefficient =T/( D) that evaluates the validity of the Stokes-Einstein law considerably varies (around a factor of 2). These results need confirmation by adding other measurements and increasing the temperature range, but on sole base of these preliminary

A C C E P T E D M A N U S C R I P T 8
experiments, it suggests that diffusion does not behave like a simple random walk like in conventional solvents. [30] Fig. 4. Variation of the peak current Ip with the square root of the scan rate (v 1/2 ) for the reduction of the ferricyanide (5 x 10 -3 mol L -1 ) in Reline) at 20°C and 60°C.

Investigations in choline-chloride/lactic acid ) eutectic (ChCl-lactic acid 1:1).

Finally, the same investigations were performed in a choline-chloride / lactic acid eutectic (1:1). As a qualitative analysis, the electrochemical signal is more defined at higher temperature especially for ferrocene oxidation (See Figure 5). Besides a possible faster electron transfer kinetics; such improvement is probably due to the concomitant decrease of the ohmic drop, a lower viscosity of the DES and a better solubility of the ferrocene in the DES at high temperature. One could notice that the peak-to-peak potential separations decrease when the temperature increases that supports a lower resistivity of the DES and/or a faster electron transfer. As obtained for the two other DESs, peak currents Ip show a linear variation with the square root of the scan rate v 1/2 that are indicative of a diffusion-controlled process in this DES as found in the two previous DES. [START_REF] Savéant | Elements of Molecular and Biomolecular Electrochemistry[END_REF] From the slopes, we derived the diffusion coefficients that are shown in Table 2. From the variation of  in Table 2, agreement with the Stokes Einstein variation is reasonable when the redox probe is ferricyanide. On the whole, measured values of the diffusion coefficients are in the same range than data reported in the literature especially if one consider that "wet" DESs are less viscous. [18][19][20][21][22][23][24][25] In Reline, diffusion was found to not follow the Stokes-Einstein law predicting that D should vary as T/. On the contrary, in Ethaline that is a less viscous DES, the Stokes-Einstein seems to prevail at least in the 20-60°C range and both for the diffusion of ferricyanide and ferrocene.

Limited or bad agreement with the Stokes-Einstein was already been reported leading to the conclusion that diffusion mechanism in DES is similar to that in the ILs [25] where a hopping mechanism or hole diffusion were proposed to explain transport inside the ILs [6][7][8]30]. The inhomogeneity of the media is often highlighted to explain the inconsistencies with the assumptions made in Stokes-Einstein equation. [30] Table 2. Diffusion Coefficients measured in the DESs. The most remarkable result of Table 2 is the low sensitivity of the diffusion coefficients to the charge of the diffusing molecule. Between ferrocene and ferricyanide, the variation of D is lower than 30% as the charge is -3 for ferricyanide and 0 ferrocene. We could also make comparisons with the self-diffusion coefficient of the choline-chloride cation (charge +1) in Ethaline that was also reported in the same order [25]. By analogy, we could remind that the diffusion of charged molecules dissolved in ionic liquids is considerably dependent on the charge. It was shown that the diffusion coefficient in ILs depends both on the delocalization of the charge on the molecule and its size. [31] For the O 2 / O 2 .-couple that is probably an extreme case, the diffusion coefficient of O 2 is 40 times larger than that of O 2 .-in common ILs. [32,33] However for larger molecules like nitrotoluene/nitrotoluene radical anion in common ILs, the neutral molecule still displays a considerable effect as diffusion coefficient of the radical anion is around 4 times smaller than the

A C C E P T E D M A N U S C R I P T

11 neutral one despite a delocalization of the charge in the aromatic ring. [31] This effect illustrates the strong ionic interactions between the redox charged species and the ions of the IL in combination with the occurrence of a possible hopping mechanism diffusion. [7] No such effects are visible in DES suggesting a much lower interaction between the solute and the solvent. Other studies are required to generalize these conclusions but from these sole results, it is clear that if DESs are comparable to ionic liquids, the transport mechanism is probably different.

Conclusion

Our experiments show that quantitative electrochemical studies could be performed in common DESs using simple cyclic voltammetry even when the redox probe is a highly charged molecule like ferrocyanide. In all considered DESs and both for ferrocene and ferrocyanide, diffusion was found as the major mass transport processes in agreement with literature results concerning similar DES.

The diffusion coefficients measured by electrochemistry were found in the same range as previously reported in previous electrochemistry investigations but also with the self-diffusion coefficients derived from NMR measurements despite different experimental times. If diffusion coefficients are smaller than those observed in classical solvents in agreement with a higher viscosity, they remain sufficiently large for use of DESs in most applications in electrochemistry.

By comparison with reports in ILs, only small variations of the diffusion coefficients were observed when changing the charge of the molecule. Other parameters as the quantity of water also need to be considered to precise how the amount of an additional component in the DES plays a role in the observed properties.

In this regard, DESs present many similarities but also differences with ionic liquids probably because of the presence of non-ionic components in the formulation of the DES. Other detailed studies are required notably with other DES to get a better view of diffusion processes in DES as mass transport is a fundamental characteristic for most applications of such solvents.

A C C E P T E D M A N U S C R I P T

  and are well adapted to studies in DESs. However, a direct extrapolation of properties in DESs from those determined in ILs is not straightforward because molecular interactions in the DESs have different origin and diversity. Although recent literature has clearly proved the potential interests of DESs, fundamental investigations in DES using electrochemistry remain limited. LeSuer et al. have reported in 2007 about the ferrocene oxidation in two choline chloride based DESs: ChCltrifluoroacetamide and ChCl-malonic acid showing the interest of the electrochemical

Fig. 1 .

 1 Fig. 1. Cyclic voltammetry of the reduction of ferricyanide (A,B) and oxidation of ferrocene (C,D) 5 x 10 -3 mol L -1 in Ethaline for different scan rates (20,50,100, 200, 500 mV s -1 ) at 20°C (A,C) and 60°C (C,D).

Fig. 2 .

 2 Fig. 2. Variations of the peak current Ip with the square root of the scan rate (v 1/2 ) for the reduction of the ferricyanide (A) and oxidation of the ferrocene (B) (5 x 10 -3 mol L -1 in Ethaline) at 20°C and 60°C.

  low solubility at 20°C and an acceptable but still low solubility at 60°C which impedes an accurate determination of the diffusion coefficient in the same experimental conditions.

Fig. 3 .

 3 Fig. 3. Cyclic voltammetry of the reduction of ferricyanide (A,B) and oxidation of ferrocene (C) 5 10 -3 mol L -1 in Reline for different scan rates (20,50,100, 200, 500 mV s -1 ) at 20°C (A) and 60°C (B,C).

Fig. 5 .

 5 Fig. 5. Cyclic voltammetry of the reduction of ferricyanide (A,B) and oxidation of ferrocene (C,D) at 5 x 10 -3 mol L - 1 in chloline chloride / lactic acid eutectic for different scan rates (20,50,100, 200, 500 mV s -1 ) at 20°C (A,C) and 60°C (C,D). Co

Fig. 6 . 10 3. 4 .

 6104 Fig. 6. Variation of the peak current Ip with the square root of the scan rate (v 1/2 ) for the reduction of the ferricyanide (A) and oxidation of the ferrocene (5 x 10 -3 mol L -1 in Choline chloride / lactic acid) at 20°C and 60°C.

Table 1 .

 1 DES used in this work.

	DES	Organic Salt	Hydrogen Donor (HBD)	Molar Ratio
	Ethaline			1:2
	Reline			1:2
	ChCl-Lactic Acid			1:1