Mathieu N Galtier

Camille Marini

Gilles Wainrib

Herbert Jaeger

Relative entropy minimizing noisy non-linear neural network to approximate stochastic processes

A method is provided for designing and training noise-driven recurrent neural networks as models of stochastic processes. The method unifies and generalizes two known separate modeling approaches, Echo State Networks (ESN) and Linear Inverse Modeling (LIM), under the common principle of relative entropy minimization. The power of the new method is demonstrated on a stochastic approximation of the El Niño phenomenon studied in climate research.

Introduction

Blackbox modeling methods for stochastic systems have a broad range of applications in physics, biology, economy or the social sciences. Generally speaking, a model of a stochastic system is a representation of the conditional distribution of the system's future given the present state (Markov models) or some part or the entire system past. There is a large variety of such stochastic predictors among which we focus on generic methods which do not depend on the type of data considered. The Auto-Regressive-Moving-Average (ARMA) models [START_REF] Box | Time series analysis: forecasting and control[END_REF] form a class of linear stochastic approximators which has led to many derivative works and is widely used in engineering applications. In particular, it covers the case of multivariate linear stochastic differential equations (SDE), or Ornstein-Uhlenbeck processes, which is the basic structure used in the Linear Inverse Modeling (LIM) theory [START_REF] Penland | Prediction of niño 3 sea surface temperatures using linear inverse modeling[END_REF]. ARMA models are generally learnt by optimizing a least squares measure of the prediction error. A notable characteristic of ARMA models is that the dimension of the underlying SDE is identical to the observable dimension of the target time series. By contrast, dynamic Bayesian networks [Murphy, 2002], with Hidden Markov Models (HMM) [Baum andPetrie, 1966, Rabiner, 1989] as the most widely employed special case, rely on hidden variables. HMM are trained by maximum likelihood schemes, typically with some version of the expectation maximization algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF], Moon, 1996]. A problem with dynamic Bayesian networks, inherited from their simpler static counterparts, is that inference (e.g. prediction) quickly becomes computationally expensive when the dependency structure of hidden variables is not particularly simple (as it is in HMMs). The Temporal Restricted Boltzmann Machine [START_REF] Sutskever | Learning multilevel distributed representations for high-dimensional sequences[END_REF], a recent addition to the spectrum of such models, is a point in case. With the advent of kernel machines in machine learning community, models based on Gaussian Processes have been designed to approximate stochastic processes [Rasmussen, 2006]. A critical point regarding these models lies in their computational complexity when working with long time series. There is also a large body of literature about online adaptive predictors, e.g.

Kalman filters [Haykin, 2005]. In this paper however we focus on non-adaptive models trained on all available training data using a batch algorithm.

Recurrent neural networks (RNNs) have also been used in various ways for approximating stochastic dynamical systems. In their basic forms [START_REF] Williams | [END_REF]Zipser, 1995, Pearlmutter, 1995], RNNs are models of deterministic dynamical systems; if trained on data sampled from stochastic sources, at exploitation time such RNNs will not propose future distributions but only a single expected mean future trajectory. RNNs represent, in principle, a promising model class because they are dense in interesting classes of target systems, implying that arbitrarily accurate models can in principle be found [START_REF] Funahashi | [END_REF]Nakamura, 1993, Sontag, 1997]. Gradient-descent based learning algorithms for RNNs are typically computationally expensive and cannot be guaranteed to converge. Since about a decade, an alternative approach to RNN design and training, now generally called reservoir computing [Jaeger andHaas, 2004, Maass et al., 2002], has overcome the problem of learning complexity. The key idea in this field is not to train all parameters of an RNN but only the weights of connections leading from the RNN "body" (called reservoir) to the output neurons. Here we will build on a particular instantiation of reservoir computing called Echo State Networks (ESNs).

Although deterministic models at the outset, neural network architectures for predicting future distributions have been variously proposed [Husmaier andTaylor, 1997, Buesing et al., 2011], or neural networks were embedded as components in hybrid models of stochastic systems [START_REF] Krogh | [END_REF]Riis, 1999, Chatzis andDemiris, 2011]. Here we propose a novel way to use RNNs in a stochastic framework based on the way stochasticity is taken into account in LIM. LIM consists in tuning both the drift and the diffusion term of an Ornstein-Uhlenbeck process to approximate a stochastic process. First, the drift is optimized to approximate the time series as if it were deterministic; then, the diffusion is chosen so that the variances of both systems are identical. LIM is widely used in climate research and stands as a simple approach giving relatively good results [Penland, 1996, Hawkins et al., 2011, Zanna, 2012, Barnston et al., 2012, Newman, 2013].

To compare two stochastic processes, and thus to define what it means to approximate a stochastic process, we use the relative entropy (also known as Kullback-Leibler divergence) [Kullback and Leibler,195 Although not a true distance, it displays many interesting properties, interpretations and relationships with other quantities such as the mutual information [START_REF] Cover | Elements of information theory[END_REF] or the rate function in large deviations theory [Ellis, 2005]. It also is computationally convenient (as opposed to the Wasserstein distance for instance), and has been widely used used in machine learning [START_REF] Ackley | A learning algorithm for boltzmann machines[END_REF], Hinton et al., 2006]. Usually, this measure is used to compare the laws of two discrete or continuous random variables, but it can also be used to compare the laws of two stochastic processes in the path space, which is at the basis of this paper. This way of measuring the difference in law between two stochastic processes amounts in performing a change of probabilities thanks to Girsanov Theorem [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF], whose applications range from mathematical finance [START_REF] Avellaneda | Calibrating volatility surfaces via relative-entropy minimization[END_REF] to simulation methods for rare events [Wainrib, 2013]. In the context of recurrent neural networks, we have already shown that the learning rule deriving from the minimization of relative entropy has interesting biological features since it combines two biologically plausible learning mechanisms [START_REF] Galtier | A biological gradient descent for prediction through a combination of stdp and homeostatic plasticity[END_REF] In this paper, we show how to train a noise-driven RNN to minimize its relative entropy with respect to a target process. The method consists two steps. First, the drift of the neural network is trained by minimizing its relative entropy with respect to the target (Section 3). Second, the noise matrix of the network is determined based on a conservation principle similarly to LIM (Section 4). We show how this approach extends the existing ESN and LIM theory in Section 5. Numerical approximations to the double well potential and to the El Niño phenomenon studied in climate research are presented in Section 6.

Model

We define here two mathematical objects that are, a priori, unrelated: a stochastic time series and an autonomous RNN made of two layers. The time series is assumed to be a sample path of an underlying stochastic process which is the modelling target. The objective is to make the RNN approximate the target process.

The target time series u is assumed to be the discretization of an n-dimensional ergodic continuous process defined on the time interval [0, T ]. The discretization step is chosen to be dt = 1 which corresponds to fixing the timescale. Imposing T ∈ N, u can be seen as a matrix in R n×T . For each t ∈ {1, .., T }, we use the notation u t for the n-dimensional vector corresponding to the value of the continuous target time series at time t. Similarly, we write δu t def = u t+1u t and δu ∈ R n×T corresponding to the previous definition (with the convention that δu T = 0).

The two-layer neural network is defined as follows. The first layer, also called retina, has n neurons, as many as the target time series dimension. We take v 0 t ∈ R n to be the activity of the retina at time t which will eventually approximate the target time series. The second layer, also called reservoir, has m neurons. Because each reservoir neuron does not directly correspond to a variable of the target, they are said to be hidden neurons. We denote the activity of the reservoir at time t by v 1 t ∈ R m . Each layer has a complete internal connectivity, recurrent connections, that is, all neurons within a layer are interconnected. The two layers are interconnected with feedforward, i.e. retina to reservoir, connections and feedback, i.e. reservoir to retina, connections, as shown in Figure 1.a. In this paper, according to a guiding principle in reservoir computing [START_REF] Lukoševičius | Reservoir computing approaches to recurrent neural network training[END_REF], the feedforward and reservoir matrices W 10 and W 11 are drawn randomly and remain unchanged. Only connections leading to retina neurons will be adapted. These are collected in W = (W 00 W 01 ) ∈ R n× (n+m) .

Retina: n neurons

Reservoir: m neurons

Reservoir activity

Retinal activity

Reservoir lateral connectivity

Retina lateral connectivity

Feedforward connectivity Feedback connectivity The activity of each layer is governed by the following differential law:

dv 0 t = -v 0 t + W 00 v 0 t + W 01 v 1 t )dt + ΣdB t dv 1 t = ǫ -lv 1 t + s(W 10 v 0 t + W 11 v 1 t ) dt (1)
where ǫ, l ∈ R + , s is a sigmoid function, e.g. tanh, that is applied elementwise, i.e. s(x) i = s(x i ), Σ ∈ R n×n is the noise matrix and B t is an n-dimensional Brownian motion.

In order to unify LIM and ESNs, we submit this architecture to certain restrictions. In particular, we choose the first layer to be linear and we choose a tanh nonlinearity in the reservoir. Later, we will make a simple choice for a numerical differentiation scheme for the same reason.

Training to minimize relative entropy

This section explains the training of the connection matrices W 00 and W 01 so that the distance between the neural network and the target time series is minimized. In other words, the drift of the neural network will be designed to match that of the target stochastic process.

Relative entropy between target and retina

We now define a quantity measuring the dynamical distance between the target time series and the retinal activity. At first sight, these two mathematical objects have a different nature: the first is a time series and the second is a dynamical system. However, we assume that there exists a dynamical system (possibly very complicated and/or with hidden variables) which has generated the target time series. Thus, we want to compute the distance between this system and the neural network. A natural measure of similarity between stochastic processes is the relative entropy or Kullback-Leibler divergence. First, we show how to compute the relative entropy between two diffusion processes with the same diffusion. Second, we apply it to computing the "distance" between the neural network and the target time series.

Relative entropy between two diffusion processes sharing the same diffusion coefficient

We now introduce the computation of the relative entropy, or Kullback-Leibler divergence, between two n-dimensional diffusion processes x and y sharing the same diffusion coefficient Σ ∈ R n×n . Let {Ω, F, P, {F t } 0≤t≤T } be a probability space, equipped with the natural filtration of the standard Brownian motion. Consider two diffusion processes (x t ) and (y t ) in R n under P , solutions of the following stochastic differential equation:

dx t = f (x t )dt + ΣdB t ( 2
)
dy t = g(y t )dt + ΣdB t (3)
where B t is a n-dimensional P -Brownian motion.

Defining and computing the relative entropy naturally follows from the application of the Girsanov theorem [Girsanov, 1960] described in chapter 3.5 of Karatzas and Shreve's textbook [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF] 1 . It provides a stochastic change of variable which makes it possible to change the drift of a diffusion process provided the underlying measure of the Brownian motion is changed accordingly. Indeed, there exists a probability measure Q and Bt a n-dimensional Q-brownian motion, such that the stochastic process (x t ) is also solution of

dx t = g(x t )dt + Σd Bt (4)
The probability measure has to be coherent with this change of drift, which is enforced through the Radon-Nikodym derivative of Q with respect to

P dQ dP = exp n i=1 T 0 Σ -1 g(x t ) -f (x t ) i dB t,i - 1 2 T 0 Σ -1 f (x t ) -g(x t ) 2 dt (5)
where Σ -1 g(x t ) -f (x t ) i and dB t,i are the i th components of Σ -1 g(x t ) -f (x t ) and dB t respectively.

For this quantity to be well-defined, the probability measure P has to be absolutely continuous with respect to Q. This is a consequence of the technical condition

T 0 E[ Σ -1 f (x t ) - g(x t ) 2 ]dt < ∞.
Given that both processes (x t ) and (y t ) are written with the same drift (3), (4), it is natural to consider the relative entropy between the processes as the relative entropy between the measures P and Q which reads

H(y|x) def = H(Q|P ) def = E Q ln dQ dP (6)
Using (5) leads to

H(y|x) = 1 2 E T 0 Σ -1 f (x t ) -g(x t ) 2 dt (7)
Strictly speaking, this quantity is not a distance since it is not symmetric. Yet it is always positive and zero only when the two drifts f and g are equal. This makes it a natural and useful measure of the similarity of two stochastic processes.

Note that although the drift g initially corresponds to the stochastic process (y t ), it is evaluated at the value x t . This is one of the main feature of relative entropy: it measures the difference of the drifts along the trajectory of one process.

It is crucial that the diffusion matrix Σ is the same in both equations ( 2) and ( 3). If it were not the case, the two measures P and Q would not be absolutely continuous and the Radon-Nikodym derivative, and a fortiori the relative entropy, would not be defined.

Application to our case

The main conceptual problem to apply the previous result is that, in practice, we do not know the continuous-time stochastic process which we assume has generated the discrete target time series. The time series is the only piece of information we have. Actually, many stochastic processes could have generated this discrete time series. We want to find the stochastic processes of the form (2), which are the more likely to have produced the time series. Let us call f the smooth function that defines the diffusion process (2) which was most likely to produce the target time series. Although we do not and will not know the explicit formula for f , we formally define the distance between the neural network (1) and the target time series as the relative entropy between the neural network and the diffusion process (2) with this particular f . We will eventually make this quantity computable.

We also need to bridge the gap between the discrete definition of the target time series and the continuous formulation of the relative entropy in ( 7). This can be done by discretizing equation ( 7) on the partition adapted to the definition of u. Assuming that the sampling of u is fine enough, we can reasonably replace the integral by a discrete sum. At this step we can also use the ergodicity property of the target time series to drop the expectation in this equation. Therefore, a first tentative of definition of the relative entropy between the target time series and the neural network is

H = 1 2T T t=0 Σ -1 f (u t ) -g(u t ) 2 + O(T -1/2 ) (8)
The term O(T -1/2 ) accounts for the ergodicity approximation when the total time window T is not infinite. However, in practice, it is not possible to have a direct access to f (u). Estimating the drift f (u) from the observation of the time-series u belongs to the class of problems called numerical differentiation. Since the seminal work of Savitzky and Golay [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF], in which the signal is first approximated by moving polynomials before differentiation, a large number of numerical methods have been introduced, from finite difference methods to regularization or algebraic methods (see [START_REF] Liu | Error analysis of jacobi derivative estimators for noisy signals[END_REF] and references therein). However, for simplicity and to rigorously relate to ESNs and LIM, we will simply approximate f (u) by a temporal difference approximation which we called δu. Recall

δu t = u t+1 -u t = t+1 t f (x s )ds + Σξ t (9)
where x t is the realization of the process from which u was sampled and the ξ t are i.i.d standard Gaussian random variables. We now assume that f is smooth enough and that the sampling of the stochastic process realization, which leads to defining the time series, is fine enough, so that the following approximation (of order 0) holds:

t+1 t f (x s )ds ≃ f (u t ) (10)
We now form the difference Σ -1 f (u t ) -g(u t ) 2 and introduce Σ -1 δu t in the computation. This leads to

Σ -1 f (u t ) -g(u t ) 2 ≃ Σ -1 f (u t ) -δu t 2 + 2 Σ -1 f (u t ) -δu t , Σ -1 δu t -g(u t ) + Σ -1 δu t -g(u t ) 2 ≃ ξ t 2 + 2 ξ t , Σ -1 f (u t ) + ξ t -2 ξ t , Σ -1 g(u t ) + Σ -1 δu t -g(u t ) 2
Let us look at what becomes each of the four lines above when taking the empirical average (or equivalently the expectation by ergodicity):

• The first line becomes 1.

• In the second line, the first term 2 ξ t , Σ -1 f (u t ) has zero mean and thus vanishes. The second term becomes 2.

• The third line is centered and thus vanishes.

• Finally, the last term is the one we want to keep in the algorithm.

To summarize, it remains

H ≃ 1 2T T t=0 Σ -1 δu t -g(u t ) 2 + O(T -1/2 ) (11)
One could also use a Taylor expansion for equation ( 10) which would add some corrective terms to (11). However, provided that the sampling of the target is fine enough, these terms can be neglected.

To be fully exhaustive, one should also look at the term O(T -1/2 ) coming from the central limit theorem for ergodic convergence, and this term may also have some contribution which depends on the connectivity, in particular from the second term in line 3: 2 ξ t , Σ -1 g(u t ) . It is zero-mean, so in the limit T → ∞ it will disappear when we take the empirical average, but its variance is in fact 4E[Σ -1 g(u) 2 ] which is not necessarily zero.

One last modification of the formal definition in equation ( 7) stems from its problematic dependence on the noise matrix Σ. Recall we intend to tune the noise to minimize the relative entropy. Given the definition (7) a trivial choice is to have an extremely strong noise: in this case the two stochastic processes are so random that they can be said to be identical. However, this is a situation which we would like to avoid: we want to find a reasonable noise matrix Σ. To correct this pathologic behavior, it seems natural to divide the formal definition (7) by the square of the operator norm of Σ -1 which we write |||Σ -1 ||| 2 (or equivalently multiplying by the square of the smallest eigenvalue of Σ). The resulting quantity simply is proportional to the original definition, but without this problem. Note that this does not mean that we restrict our approach to unit norm noises matrices as is shown later.

This finally leads to a definition corresponding to equation ( 11) without the various error terms. Besides, we also want to take into account a regularization term. This leads to defining an analog to regularized relative entropy between the target time series u and our neural network (1) (respectively corresponding to x and y in (2) and ( 3)) as

H u (W) def = 1 2T T t=0 S -1 -u t + W 00 u t + W 01 u 1 t -δu t 2 + α 2 2 S -1 W 2 (12)
where

S -1 = Σ -1 |||Σ -1 ||| and u 1 ∈ R m×T is the solution of δu 1 t = ǫ -lu 1 t + s(W 10 u t + W 11 u 1 t )
governing the reservoir when it is fed by the target time series (as opposed to the retinal activity). The definition of δu 1 ∈ R m×T is similar to that of δu. In the following, we abusively call H u (W) the relative entropy between the target time series and the network.

The second term in ( 12) is a regularization term which ensures decent generalization properties of the following algorithm. It does not change the qualitative meaning of the notion of relative entropy but penalizes the networks with high connection strength. The reason why it contains a multiplication by Σ -1 will become clear later. In short, it will make the system decoupled and we will be able to compute subsequently the connections first and the noise matrix second.

One can also understand the definition (12) without referring to the analogy with relative entropy, in the special case Σ = I d and α = 0. Indeed, as shown in Figure 2, it corresponds to integrating the distance between the derivative of the target time series and the flow of the retinal activity along the trajectories of the target time series. Therefore, it is obvious that if this quantity is null then the vector field of the retinal activity will be tangent to the target time series derivative. Therefore, initializing the network with the value of the target time series at time t = 0, will lead the network to reproduce precisely the target time series. Observe that the relative entropy (12) (even more the rigorous definition ( 7)) are quite similar to the quantity usually minimized for prediction [START_REF] Williams | Gradient-based learning algorithms for recurrent networks and their computational complexity[END_REF], Pearlmutter, 1995, Bishop, 2006]. The main difference is that this is the integral of the distance between the derivative of the target time series and activity variable, instead of the mere distance between target and activity.

Note that the relative entropy (12) can easily be shown to be proportional to the negative of the log-likelihood of the target given the neural network. This close relationship shows that, in this case where the stochastic processes are diffusion processes, minimizing the relative entropy is rigorously equivalent to maximizing the log-likelihood.

Gradient of the relative entropy

Now, the idea is to compute the gradient of H u (W) with respect to retinal and feedback connectivities. This will be useful because the gradient cancels out at the minimum of the relative entropy. This gives us a useful way to compute the connectivity W * = (W * 00 W * 01 ) ∈ R n×(n+m) such that the network best approximates the target stochastic process. Because the relative entropy ( 12) is quadratic in W, it is convex. Thus, it has a single critical point which is a global minimum.

For readability, we introduce the continuous nn-dimensional function defined on [0, T ] by

ζ(W) = W 00 u + W 01 u 1 -δu -u (13)
such that

H u (W) = 1 2 T t=0 S -1 ζ t (W) 2 + α 2 2 S -1 W 2 (14)
The gradient will be computed in an operator framework, i.e. as the differential of the relative entropy. Observe that the differential is a linear operator such that

d W H u (J) = 1 2 d W T t=0 S -1 ζ t , S -1 ζ t + α 2 S -1 W, S -1 W (J) = T t=0 1 2 d W S -1 ζ t , S -1 ζ t (J) + α 2 S -1 W, S -1 J = T t=0 S -1 ζ t (W), S -1 d W ζ t (J) + α 2 S -1 W, S -1 J = T t=0 SS ′ ) -1 ζ t (W), d W ζ t (J) + α 2 (SS ′ ) -1 W, J (15)
Because ζ(W) is affine in W 00 and W 01 , it appears that

d W 00 ζ(J) = Ju and d W 01 ζ(J) = Ju 1 (16)
In both cases, we have a differential of the form

d W ζ(J) = Jb, with b = u (resp. b = u 1 ) is a matrix of R n×T (resp. R m×T ). Observe that ∇ W H uij = d W H u (E ij
) where E ij is the canonical matrix made of zeros except at position ij where it is one. Then,

d W H u (E ij ) = T t=0 (SS ′ ) -1 ζ t (W), E ij b t + α 2 (SS ′ ) -1 W, E ij = T t=0 (SS ′ ) -1 ζ t (W) i b jt + α 2 {(SS ′ ) -1 W} ij = (SS ′ ) -1 ζ(W)b ′ ij + α 2 {(SS ′ ) -1 W} ij which leads to ∇ W H u = (SS ′ ) -1 ζ(W)b ′ + α 2 (SS ′ ) -1 W
Using the definition of ζ leads to

(SS ′ )∇ W H u = -(δu+u)u ′ (δu+u)u 1 ′ +(W 00 W 01 ) α 2 I d + uu ′ uu 1 ′ u 1 u ′ u 1 u 1 ′ (17)
The global minimum can be computed as ∇ W * H u = 0. Since we have assumed that S is a full rank matrix, this leads to the following formula

W * = (δu + u)u ′ (δu + u)u 1 ′ α 2 I d + uu ′ uu 1 ′ u 1 u ′ u 1 u 1 ′ -1 = (δu + u) u u 1 ′ α 2 I d + u u 1 u u 1 ′ -1 (18)
It is interesting to observe that the solution does not depend on the noise matrix Σ. It is this property which makes it possible for the problem to be decoupled in two parts: (i) computation of W * and (ii) computation of Σ.

4 Computing the noise with a conservation principle

This section is devoted to computing the noise matrix Σ such that the neural network matches the statistics of the target time series. Obviously the choice of a simple additive noise (i.e. Σ does not depend on v nor t) in system (1), restricts the class of target system the neural network can approximate accurately. Although the match will not be perfect, we will see the method provides a reasonable and (more crucially) coherent noisy neural network approximating the statistics of the target.

At first sight, it may seem appropriate to choose Σ so that the covariance of the retinal activity matches that of the target. However, the non-linearity in the reservoir makes it seemingly impossible to compute analytically the covariance of system (1). Therefore, we have not been able to use this idea to fix Σ.

Another method consists in using a generalized fluctuation dissipation relation. Penland and Matrosova [START_REF] Penland | A balance condition for stochastic numerical models with application to the el nino-southern oscillation[END_REF] have detailed a method to use a conservation principle to link the correlation of the activity, the flow of the retina and the matrix ΣΣ ′ . Following the lines of their derivation, we generalize their approach to the case of reservoirs.

The generalized fluctuation dissipation relation is based on the Fokker-Planck equation [Risken, 1996] of the neural network (1). Any stochastic differential system can be described equivalently by a sample path dynamics governed by (1) or a Fokker-Planck equation which governs the evolution of the probability density function. It corresponds to the Eulerian description of the original stochastic differential equation. It can intuitively be understood as a balance of how much goes in and out of a small box centered on v, taking into account both a drift and a diffusion mechanism. In our case, it takes the form of the following partial differential equation.

∂p(v 0 , v 1 , t) ∂t = -div -v 0 + W 00 v 0 + W 01 v 1 ǫ -lv 1 + s(W 10 v 0 + W 11 v 1 ) p(v 0 , v 1 , t) + 1 2 ∆ ΣΣ ′ 0 0 0 p(v 0 , v 1 , t) (19)
where div is the divergence operator, i.e. div(x) = i ∂x i ∂v i which corresponds to the drift, and ∆ is the Laplacian operator, i.e. ∆J = i,j ∂ 2 J ij ∂v i ∂v j , which corresponds to the diffusion. Note that this Fokker-Planck equation is independent of the underlying choice between Itô or Stratonovich noise in the initial system (1), because we assumed Σ does not depend on the activity of the network.

A conservation principle about the moments of the stochastic process v can easily be derived from equation (19). Indeed, multiply it by {v 0 t } p {v 0 t } q (the components p and q of the activity v 0 t ) and integrate over the entire domain. We can use the integration by part formula several times to get

∂E[{v 0 t } p {v 0 t } q ] ∂t = E -{v 0 t } p + {W 00 v 0 t } p + {W 01 v 1 t } p {v 0 t } q + E -{v 0 t } q + {W 00 v 0 t } q + {W 01 v 1 t } q {v 0 t } p + ΣΣ ′ (20)
Given that the target is ergodic, it is natural to assume that its approximation also has this property. Therefore, it is legitimate to replace the expectations by time integrals over [0, T ] divided by T in the previous equation. In a matrix formalism, this reads

v 0 T v 0 T ′ -v 0 0 v 0 0 ′ T = -v 0 + W 00 v 0 + W 01 v 1 v 0 ′ + v 0 -v 0 + W 00 v 0 + W 01 v 1 ′ + ΣΣ ′
A careful inspection of the terms shows that the left hand side is negligible when T is large enough (which will always be the case in practice). Thus, we will drop this term for simplicity (although it would not pose any problem to take it into account). Notice the correlations terms v 0 v 0 ′ and v 0 v 1 ′ in the previous equation. Recall our initial wish to choose Σ so that both neural network and target second order moments are matched. Although we could not directly implement this wish, we are now able to replace the network correlation terms by the observed moments of the target in the present formulation. This ansatz leads to the following generalized fluctuation dissipation relation [START_REF] Penland | A balance condition for stochastic numerical models with application to the el nino-southern oscillation[END_REF]:

ΣΣ ′ = I d -W 00 uu ′ T -W 01 u 1 u ′ T + uu ′ T I d -W 00 - uu 1 ′ T W 01 (21)
This equation can be seen as a coherency requirement between drift and diffusion of the network and second order moments of the target. Fortunately, the derivation of the drift leads to an explicit equation ( 18) independent of the matrix Σ. Therefore, the previous equation can be used to characterize Σ based on the knowledge of W. Given that the square root of matrix is not injective, there are several choices for the matrix Σ. They all correspond to an ambiguity on the sign of its eigenvalues. We arbitrarily pick one of them and have thus found a coherent noise matrix.

Comparison with existing methods

Our approach takes selected features of two existing methods for approximation, ESN and LIM, and unifies them in a mathematical framework. This section is devoted to clarifying the links with these two methods. 2): the inputs and outputs are different layers. However, when the network is set to the prediction mode the value of the inputs is copied from the outputs, which corresponds to closing up the loop as in Fig. 1. b) Architecture of LIM (see section 5.1): there is no reservoir and only a retina.

Adding nonlinearities to Linear Inverse Modeling

The proposed method rigorously extends LIM by adding non-linearities to the model dynamics.

Linear inverse modeling consists in designing a multidimensional linear stochastic (also called Ornstein-Uhlenbeck) process which reproduces a target multivariate time series [Penland andSardeshmukh, 1995, Penland, 1996]. Naturally, the dimension of the approximating process is identical to the dimension of the target. A pervasive idea in machine learning is to consider additional (often called hidden) variables which will help the reconstruction of the target multivariate time series. In the framework presented in this paper, they correspond to the neurons in the reservoir, while the linear dynamical system analogous to LIM corresponds to the neurons in the retina. Due to the nonlinearity of the dynamics of these additional variables or neurons, the present framework is as a non-linear extension of LIM. Actually, it turns out to be a surprisingly simple extension since the same formula can be used to calculate the linear matrix in LIM and the rectangular matrix which combines the additional variables to improve the retina's predictions.

More precisely, LIM consists in finding the matrices B, Q ∈ R n×n such that the following dynamical system reproduces the target time series u.

v0 t = Bv 0 t + Q Ḃt
where Ḃt is a white noise and √ is the matrix square root (i.e. Q is the covariance of √ Q Ḃt ). Note that we have intentionally used the same variable v 0 for the LIM and the activity in the retina.

Based on the explicit expression of the Green function of a linear system [Risken, 1996], B can be caracterized as follow

B = 1 τ ln(u .+τ u ′ (uu ′ ) -1 ) ( 22
)
where ln is the logarithm for matrices and u .+τ is the matrixes whose column number t is u t+τ and τ is an integer usually equal to 1 (depending on the intrinsic timescales of the target). Assuming τ = 1 in the following, B can be written:

B = ln(I d + (u .+1 -u)u ′ (uu ′ ) -1 ) = +∞ k=1 (-1) k+1 k (u .+1 -u)u ′ (uu ′ ) -1 k
Assuming appropriate sampling of u such that δu ≪ 1, we can reasonably truncate B at first order. This leads to the following approximation: B ≃ δuu ′ (uu ′ ) -1 .

To see the link with ESNsto, we must consider the case m = 0 in eq.1, so that W = I d + B. From the above approximation, it follows 18) with m = 0 and α = 0. Given the definition of both systems, it is clear that the two methods are identical in this restricted case. However, when the sampling is not fine enough, equation ( 22) may lead to solutions than differ from the log-free equation ( 18). One then could also imagine a variation of ESNs using a matrix log when dealing with badly sampled data, but this is beyond the scope of this paper.

W ≃ I d + δuu ′ (uu ′ ) -1 ≃ (δu + u)u ′ (uu ′ ) -1 , which corresponds to W * in (
Concerning the treatment of noise, we can observe that it is strictly the same with or without additional reservoir neurons.

Adding noise to Echo States Networks

The proposed algorithm rigorously generalizes ESNs to a stochastic framework. Indeed, the formula for the connectivity in equation ( 18) is identical to the solution given by applying the classical deterministic ESN method [Jaeger andHaas, 2004, Lukoševičius andJaeger, 2009].

To explain further the equivalence between this formalism and classical ESNs, we now introduce the ESN formalism in its original form as summarized in Figure 3.a. Let us setup an ESN of reservoir size m as a one time step predictor of the input u, where the input and the output have, naturally, the same dimension n, and the model includes direct connections from the input to the output. In fact, the initial setup of echo state network makes a distinction between input and output, see Figure 3.a; whereas the model introduced in section 2 only has a retina, see Fig. 1. However, ESNs can also be run in a generative mode, where the current output becomes the next input, closing the loop between the two. This closed loop system is precisely the same as the two layer network of Section 2, where the joined input/output nodes become the retina with activations v 0 , and the reservoir remains with activations v 1 . As summarized in Figure 3.a, connections from the input to the reservoir correspond to W 10 , internal reservoir connections to W 11 , output connections from the reservoir to W 01 , and connections from input directly to the output after closing the loop become the recurrent connections in the retina W 00 .

Classical ESNs are discrete-time systems, as opposed to our continuous-time approach. Yet the two methods are closely linked: ESNs correspond to a time-discretized version of (1).

v 0 t+1 = W 00 v 0 t + W 01 v 1 t v 1 t+1 = s(W 10 v 0 t + W 11 v 1 t ) , (23)
where ( 1) is discretized using Euler's approximation and the discretization step is taken to be equal to 1 and ǫ = l = 1.

Training such a setup to minimize a squared error on the input prediction T t=0 u t+1v 0 t 2 precisely corresponds to learning connections W 00 and W 01 according to the equation ( 18). Indeed, observe that {δu+u} t = u t+1 is the prediction of the input (which corresponds to the target signal in this setup). Taking u as input and u 1 as teacher-forced reservoir activations, the equation ( 18) turns out to be the ridge regression equation generally used in [START_REF] Lukoševičius | Reservoir computing approaches to recurrent neural network training[END_REF].

The treatment of noise proposed in this paper is a new contribution to the ESN theory. In that sense, this paper consists in an extension of ESNs as time series approximators to stochastic ESNs as stochastic process approximators.

Numerical simulations

This section shows two examples of application of the proposed ESNsto algorithm described in algorithm 1.

Algorithm 1 ESNsto learning

# Initialization: Components of W 01 and W 11 are drawn randomly following a normal law. The two matrices are then rescaled get desired spectral radii.

u 1 0 , v 0 0 , v 1 0 ← 0 # Collect reservoir states: while t < T : do u 1 t+1 = (1 -ǫl)u 1 t + ǫl s(W 10 u t + W 11 u 1 t ) end while # Learn connections: (W 00 , W 01 ) ← (δu + u) u u 1 ′ α 2 I d + u u 1 u u 1 ′ -1
# Compute noise connections:

Σ = I d -W 00 uu ′ T -W 01 u 1 u ′ T + uu ′ T I d -W 00 -uu 1 ′ T W 01 # Simulate ESNsto: while True do v 0 t+1 = W 00 v 0 t + W 01 v 1 t + Σ randn(n) v 1 t+1 = (1 -ǫlv 1 t ) + ǫs(W 10 v 0 t + W 11 v 1 t ) end while
It is important to realize that the goal here is not to approximate or predict a time series, but rather to approximate a stochastic process. Because a single stochastic process can have different realizations, an approximation of such a process should not aim at reproducing a given path. In this sense, we are not dealing with classical prediction tasks and we should exlusively focus on building a system that reproduces the law of the target stochastic process. As a consequence, we can not compare ESNsto with classical ESN since they do not approximate the same mathematical objects.

We are going to compare the performances of LIM with that of ESNsto. LIM belongs to several of the different classes of approximators that we mentioned in the introduction: it is a Gaussian process, a multivariate autoregressive process of order one and an Ornstein-Uhlenbeck process. Besides it can be easily compared to ESNsto, since the latter generalizes the former and LIM simply is an ESNsto with 0 neurons in the reservoir. Establishing a complete benchmark of the different methods for stochastic processes approximation is beyond the scope of this paper. However, we point out the low computational complexity of learning, which is independent of the length of the time series and mainly governed by the inversion of a positive semi definite square matrix of size n. This is lower than the complexity of Hidden Markov Models or Gaussian Processes for long time series.

To compute the relative entropy we have used a classical cross-validation framework. This means that we have divided the target time series in k blocks. Then for each block, we have computed the connectivity and noise matrices on the remaining k -1 blocks, and evaluated the value of the relative entropy on the selected block. This provides a robust way to prevent overfitting. This numerical section is only a proof of concept. The (hyper)parameters of the networks (such as the spectral radius of W 11 ) have been coarsely tuned, although they could significantly improve the approximations if they were set carefully. The main reason for our negligence is that we want to show that an off-the-shelf ESNsto model is better than LIM and does not require deep knowledge or large effort of neural networks tuning.

The first example we consider is a widely considered toy model: the target is generated by a noisy particle living in a double well. The second example, devoted to climate modeling, will show how to approximate the El Niño phenomenon in the tropical Pacific ocean.

The noisy double-well

The double well example explores a basic form of non-linearity. It illustrates the significant improvement brought by reservoir neurons in dealing with non-linearities. We consider a synthetic example where the data are generated as the solution of the stochastic differential equation corresponding to a particle in an energy landscape made of two different wells, as shown in Fig. 4(a). More precisely, the target is a one-dimensional process described by

du = -∇ u Edt + σdB t ( 24
)
where

∇ u E = 1 if -1 < u < 0 or u > 1 -1 else
is the gradient of the function described in Fig. 4(a) at point u. The typical behavior of such a system is illustrated in Fig. 4(b). Roughly speaking, the particle jumps from one well to the other after random durations. Informally, each well can be said to be an attractor. We now compare the approximations of this time series based on LIM and the ESNsto. We see in Fig. 5(a) that increasing the number of neurons in the reservoir improves the relative entropy defined by (12). With no neurons in the reservoir, which corresponds to LIM, the relative entropy is approximately 2.9. After a sharp decrease for the first dozens of additional neurons, the relative entropy slowly decreases when the number of neurons increases to finally reach a value close to 2.65 for m = 150 neurons. Note that the simple numerical differentiation method that we have used imposes a lower bound on the relative entropy shown here. Indeed, it is easy to observe that the relative entropy between system (24) and itself using definition (12) is 2.45. We believe that the gap between the ESNsto and the optimal value will decrease with additional neurons in the network, but may not vanish due to the over-fitting issue mentioned later. This suggests that appropriately choosing the regularization parameter is crucial for optimal accuracy. We can also observe that the variance of the relative entropy value, corresponding to different random realizations for the connections in the reservoir W 11 and W 10 , decreases with the number of neurons: due to better averaging, large reservoirs are less dependent on the realization defining their weights.

Running the LIM and ESNsto networks post-learning shows different qualitative behaviors as displayed in Fig. 5 W 11 is drawn according to a normal law and rescaled so that W 11 = 1 (operator norm), W 10 is drawn uniformly in [-1.1], α = 1. The network simulations and learning were done during T = 10 5 time steps. Note that there is a washout time interval of 500 at the beginning of the learning (and reconstruction) in order to remove transient effects.

Significant qualitative improvements made possible by having neurons in the reservoir can also be seen by empirically measuring some statistical quantities of the target, LIM and ESNsto runs, as shown in Fig. 6. The first three figures 6(a), 6(b) and 6(c), show that LIM is failing to reproduce the bimodal distribution of the target corresponding to the two attractors: LIM has a unimodal Gaussian-like distribution, whereas ESNsto is able to reproduce the bimodality thanks to the non-linearities in the reservoir.

In Fig. 6(d), it is shown that the distribution of the times spent in each attractors between two jumps is irrelevant for LIM whereas it is similar between target and ESNsto. This is also reflected in the transition rates which is approximately 0.0019 for the data, 0.0016 for the ESNsto and 0.0075 for the LIM. However, it is to be noticed that an increase in the number of neurons beyond 150 neurons leads to a decrease of the transition rate (not shown). This underlines a drawback of the method for non-ergodic time series which exhibit significant noisy flucutuations: the network tries to put in the connectivity as much variability as possible; the noise term is simply taking care of the left-overs. Therefore, the noise induced transitions in the target are not only modeled by the diffusion term in the neural network but also by the drift. This effect will vanish if the learning time series has enough ergodicity so that the noisy behavior is averaged out when computing the drift. When there is a limited amount of time steps available, a better numerical differentiation scheme may improve the approximation accuracy since it would filter out noise before asking the drift to approximate it.

El Ni ño phenomenon

In this section, we focus on approximating the geophysical process El Niño. It corresponds to a large warming of Sea Surface Temperature (SST) in the eastern equatorial Pacific, occurring irregularly every 2 to 7 years, and having a broad impact of the global climate [Trenberth, 1997, Deser et al., 2010]. As many other geophysical processes, its dynamics evolves on different interacting timescales. It is in particular strongly linked to atmospheric processes evolving at shorter time scales with a nonlinear behaviour. Penland showed that the evolution of 3-month running mean SST anomalies in the tropical Indo-Pacific, and thus of El Niño, are well approximated by a LIM, where the rapidly varying nonlinear processes are parameterized as a stochastic forcing of the slower system [Penland, 1996]. One commonly used index of the El Niño phenomenon is the Niño 3.4 index (N34 index), defined as the averaged of SST anomalies between 5°S-5°N and 170°W-120°W. Using the definition of [Trenberth, 1997], an El Niño event is said to occur when the N34 index, smoothed with a 5-month running mean, exceed 0.4°C for 6 months or more. Fig. 7 shows the N34 index (top) and the regression of SST anomalies onto this index (bottom), indicating the warming in the eastern equatorial Pacific associated with a positive N34 index.

The target time series considered here are the N34 index, smoothed with a 3-month running mean, and the 10 first principal components (PCs) of the empirical orthogonal function of 3month running mean SST anomalies in the IndoPacific region (30°S-30°N, 40°E-70°W). These 10 PCs represent 80% of the total variance of monthly IndoPacific SST anomalies, and they are used instead of considering directly the SST anomalies at each grid point in the IndoPacific to reduce the dimensionality of the system. The data come from the HadISST1 SST dataset [START_REF] Rayner | Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century[END_REF], constructed from in situ SST observations and satellite derived estimates and available from 1870. We used data from 1870 to 2011. Our target times series contain thus 1704 time steps (corresponding to 1704 months or 142 years) and a washout period of 240 time steps is used at the beginning of the learning to remove transient effects. Here, we compare the approximations of the N34 index based on LIM and ESNsto. A crucial parameter for the success of ESNsto was the choice of ǫ = 0.1 in equation ( 1). It is known that this parameter controls the speed of the reservoir [START_REF] Jaeger | Optimization and applications of echo state networks with leaky-integrator neurons[END_REF]. In our case, the reservoir needed to evolve, not at the scale of months (which would have corresponded to ǫ = 1), but rather at longer time scales to be helpful in reconstructing the dynamics. Fig. 8 shows the relative entropy of the system as a function of the number of neurons and the ridge regularization parameter α. Without any regularization, the relative entropy increases with the number of neurons. However, for α > 10, the relative entropy decreases with the number of neurons and with α, suggesting that ESNsto leads to a better approximation than LIM when using strong regularization. This can be interpreted as overfitting in the case of weak regularization. Adding regularization penalizes the accuracy on the training dataset (not shown) but significantly improve the generalization on the test dataset, as observed in Fig. 8 This choice of m and α is motivated by Fig. 8. The system is simulated forward for 10 5 time steps, corresponding to more than 8000 years.

As shown in Fig. 9(a), the spectrum of the N34 index is closer to the target when approximated by ESNsto than by LIM. At time scales longer than 4-5 years, the latter shows too much variability. The spectrum obtained with ESNsto also shows a too high variability, but less than the LIM and only at time scales longer than 7 years.

The distributions of the N34 index based on the LIM and ESNsto simulations are compared with the targeted distribution in Fig. 9(b). Due to the low number of target samples, it is hard to determine from the figure which distribution is closer to the target. A Kolmogorov-Smirnov test is used to determine whether the simulated distributions differ from the targeted one. The p-values in the case of LIM and ESNsto are respectively 0.45 and 0.62, meaning that, in both cases, we cannot reject the null hypothesis that the simulated and targetted N34 indices are drawn from the same distribution. The larger p-value in the case of ESNsto indicates a stronger evidence against the null hypothesis.

Fig. 9(c) shows the distributions of the time interval between two El Niño events. Again, a Kolmogorov-Smirnov test is used to estimate if the simulated distributions differ from the targeted one, and gives a p-value of 0.76 for LIM and 0.98 for ESNsto. The p-value almost equals 1 in the case of ESNsto, suggesting a very good accuracy of our model to reproduce the some aspects of the dynamics governing El Niño events.

Discussion

We have shown how to design a recurrent neural network to reproduce a target stochastic process. By doing so, we have introduced a rigorous mathematical derivation which unifies ESNs and LIM under the general principle of relative entropy minimization. Finally, we have shown how the proposed system outperforms LIM, even when the parameters are only coarsely tuned, on a simple synthetic task and on a climate example of well-known importance.

We have observed that the system is prone to over-fitting, which forced us to use large regularization parameters. Indeed, the sequential computation of connectivity matrix followed by noise matrix, implies that noise only takes care of left-overs. The connectivity matrix will try to encode as much of the signal as possible, even some part of the inherent noise. This is problematic in applications where the target process is significantly noise driven. However, when the number of time steps of the target time series is large enough to have a good ergodic approximation or if we improve the numerical differentiation scheme (e.g. using the Savitzky-Golay algorithm), we believe this drawback will vanish.

Figure 1 :

 1 Figure 1: Structure and main notations of the neural network described in Section 2.

Figure 2 :

 2 Figure 2: Illustration of the computation of the relative entropy in the case Σ = I d . It corresponds to integrating all the green bars along the trajectory of the target time series. The solid blue line corresponds to the target time series. The red arrows correspond to the speed vectors of the target time series. The black arrows correspond to the vector field of the retinal neural network.

Figure 3 :

 3 Figure3: Architecture of related existing models. To be compared with Fig.1. a) Architecture of an ESN (see section 5.2): the inputs and outputs are different layers. However, when the network is set to the prediction mode the value of the inputs is copied from the outputs, which corresponds to closing up the loop as in Fig.1. b) Architecture of LIM (see section 5.1): there is no reservoir and only a retina.

Figure 4 :

 4 Figure 4: Illustrations of the double well process (24). (left) Energy landscape. (right) Trajectory of the process. The simulation was done using Euler-Maruyam algorithm with σ = 0.7, dt = 0.1 during 10000 time steps. These parameters are kept constant for the next simulations.

  (b) and 5(c). As opposed to LIM, the ESNsto reproduces patterns of noiseinduced jumps between two attractors. number of neurons (a) Relative entropy vs number of neurons. time-steps (b) LIM trajectory. time-steps (c) ESNsto trajectory.

Figure 5 :

 5 Figure 5: (a) Relative entropy as a function of the number of neurons in the reservoir. For any number of neurons, we launched 10 simulations, to take into account the variability of reservoir connections, and measured the relative entropy according to equation (12). The top (resp. bottom) curve is the maximum (resp. minimum) value of the entropy among the different trials. The middle curve is the mean. The relative entropy roughly decreases from 2.9 to 2.65 which is enough to qualitatively change the behavior of the neural networks as shown in the left and middle picture. (b) Result of the LIM simulation (i.e. m = 0) learned from the data in Fig.4(b). (c) Result of the ESNsto simulation with m = 100 neurons. Parameters used: ǫ = 1, l = 1, s(x) = tanh(2x),W 11 is drawn according to a normal law and rescaled so that W 11 = 1 (operator norm), W 10 is drawn uniformly in [-1.1], α = 1. The network simulations and learning were done during T = 10 5 time steps. Note that there is a washout time interval of 500 at the beginning of the learning (and reconstruction) in order to remove transient effects.

  State distribution for the data.

  State distribution for the LIM.

  Distributions of the durations between jumps.

Figure 6 :

 6 Figure 6: Distribution of the states of the data (a), LIM (b) and ESNsto (c). (d): Compared histograms of the duration spent in a well between two jumps. The parameters used are the same as for Fig.5.

Figure 7 :

 7 Figure 7: Niño 3.4 index, smoothed with a 3-month running mean (top) and associated SST anomalies (bottom).

  (b).

Figure 8 :

 8 Figure 8: Relative Entropy as a function of the number of neurons m and the ridge regularisation parameter α. Parameters used: ǫ = 0.1, l = 1/2, s(x) = tanh(2x), W 11 is drawn according to a normal law and rescaled so that W 11 = 1 (operator norm), W 10 is drawn uniformly in [-1.1], α = 1. Washout = 240. The cross-validation was done with k = 10 blocks. Each value of the relative entropy corresponds to the average over 10 realizations of the weight matrices. We now compare the simulations based on LIM and ESNsto with m = 500 and α = 3000.

  frequency in cycles per month (a) Spectrum of the N34 index. Frequency in cycle per month. (b) Distribution of the N34 index (c) Distribution of the time interval between two El Niño events

Figure 9 :

 9 Figure 9: Spectrum (left) and distribution (middle) of the N34 index, smoothed with a 3-month running mean. Distribution of the time interval between two El Niño events (right). Black corresponds to the targeted time series, blue to LIM, and red to ESNsto. Parameters are identical to those of Fig. (8).

we apply the theorem 5.1 with, in their notations, Xt = Σ -1 g(xt) -f (xt) , Wt = Bt, Wt = Bt and ZT is the Radon-Nikodym derivative of Q with respect to P according to equation 5.4.

Possible extensions of this theory could include a proper treatment of the case of badly sampled data as well as the generalization of the method to space dependent diffusion coefficients. Finally, an important step in increasing accuracy of these networks would be to identify an appropriate automatic tuning of the hyper parameters of the network (e.g. the spectral radius of the reservoir connections).

Acknowledgment

The authors would like to thank Mantas Lukosevicius for helpful discussions. MNG was funded by the Amarsi European Project.


Chiara Sinito

Mark J Fernée

Serguei Goupalov

Philippe Tamarat

Brahim Lounis

Spectroscopic signature of shape anisotropy in single CdSe nanocrystals

Spectroscopic signature of shape anisotropy in single CdSe nanocrystals

C. Sinito 1 , M. J. Fernée 1 , P. Tamarat 1 and B. Lounis 1 1 Université de Bordeaux, LP2N F-33405 Talence, France E-mail: chiara.sinito@u-bordeaux1.fr

The understanding of the photophysical properties of semiconductor nanocrystals is of utmost importance for the development of their applications in the fields of laser technology, biological labeling, nanoelectronics and quantum information. Single nanocrystal spectroscopy and the theoretical model developed by Efros, Rosen et al. [1] have provided a deeper understanding of the properties of the band edge exciton fine structure in spherical CdSe NCs with ellipsoidal perturbations. We study single zinc blende CdSe nanocrystals for the high symmetry of their cubic structure, which makes them very sensitive to deformations around the spherical shape [2]. We have developed a photoluminescence excitation technique (PLE) to probe the whole band edge exciton fine structure. We observe both split 1 L and 1 U states, which we attribute to shape anisotropy. The polarizations of the doublets are anticorrelated as predicted by the theory [3]. The split states are not orthogonally polarized, which can be interpreted as the result of non-orthogonal perturbations of the spherical shape. Moreover, states coupling induced by an external magnetic field gives oscillator strength to forbidden transitions, allowing to resolve for the first time all the eight states of the band edge exciton fine structure, as shown in Fig. 1 for an elongated nanocrystal. Optical spectroscopy therefore produces the full spectral fingerprint of single nanocrystals, which is essential to test long standing theoretical models.

Fig. 1

 1 Fig. 1 Left : band edge exciton fine structure splitting as a function of the aspect ratio for zinc blende CdSe nanocrystals (dotted lines correspond to dark states, solid lines indicate bright states). Right : PLE scan at 2 Kelvin and 6 Tesla showing the whole fine structure of the band edge exciton for a single elongated nanocrystal.


Vesna Andova
email: vesna.andova@gmail.com

František Kardoš
email: frantisek.kardos@labri.fr

Riste Škrekovski
email: skrekovski@gmail.com

Fullerene Graphs and Some Relevant Graph Invariants

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Fullerenes are polyhedral molecules made entirely of carbon atoms. They come in wide variety of sizes and shapes. The most symmetric is the famous buckminsterfullerene, C 60 , whose discovery in 1985 marked the birth of fullerene chemistry [40].

The name was a homage to Richard Buckminster Fuller, whose geodetic dome it resembles. Since 1991, when the Science magazine pronounced the buckminsterfullerene for the "Molecule of the year", fullerenes have been attracting attention of diverse research communities. The experimental work was paralleled by theoretical investigations, applying the methods of graph theory to the mathematical models of fullerene molecules called fullerene graphs.One of the main driving forces behind that work has been a desire to identify structural properties characteristical for stable and reactive fullerenes, i.e., for fullerene isomers verified in macroscopic quantities. Ante Graovac also found fullerenes as an interesting topic to study [START_REF] Ghorbani | An upper bound for energy of matrices associated with an infnite class of fullerenes[END_REF][START_REF] Kaufman | Graph drawing algorithms versus molecular mechanics geometries in fullerenes[END_REF][START_REF] Plestenjak | Generating fullerenes at random[END_REF], he connected the well known topological indices with the fullerene structures [START_REF] Fowler | The Estrada index and fullerene isomerism[END_REF][START_REF] Graovac | Distance property of fullerenes[END_REF].

Fullerene graphs

Fullerenes can also be seen as graphs, vertices represent atoms, and edges represent bonds between atoms. A fullerene graph is a 3-connected 3-regular planar graph with only pentagonal and hexagonal faces. By Euler's formula, it follows that the number of pentagonal faces is always twelve. Grűnbaum and Motzkin [START_REF] Grünbaum | The number of hexagons and the simplicity of geodesics on certain polyhedra[END_REF] showed that fullerene graphs with n vertices exist for all even n ≥ 24 and for n = 20. Although the number of pentagonal faces is negligible compared to the number of hexagonal faces, their layout is crucial for the shape of a fullerene graph. There are fullerene graphs where no two pentagons are adjacent, i.e., each pentagon is surrounded by five hexagons. Those fullerene graphs satisfy the isolated pentagon rule or shortly IPR, and they are the most stable fullerene compounds [START_REF] Kroto | The stability of the fullerenes C n[END_REF]. If all pentagonal faces are equally distributed, we obtain fullerene graphs of icosahedral symmetry, whose smallest representative is the dodecahedron. The dodecahedron is the only icosahedral fullerene that does not satisfy the IPR. On the other hand, if the pentagonal faces are grouped in two clusters by six, we obtain nanotubical fullerene graphs.

Icosahedral fullerene graphs

The common feature of all icosahedral fullerenes is their geometrical shape. The simplest icosahedral fullerene graph is the dodecahedron, C 20 . Caspar and Klug [START_REF] Caspar | Viruses, Nucleic Acids and Cancer[END_REF] and Coxeter [START_REF] Coxeter | Virus macromolecules and geodesic domes[END_REF] suggested a method that works with the icosahedral fullerene graph dual: geodesic domes, i.e., triangulations of the sphere with vertices of degree 5 and 6.

Goldberg [START_REF] Goldberg | A class of multi-symmetric polyhedra[END_REF] observed that all icosahedral fullerene graphes can be obtained by mapping (a part of) the hexagonal grid onto the triangular faces of an icosahedron.

He also showed that the number of vertices n in a polyhedron of icosahedral symmetry can be determined by two integers i and j by the following equation, conveniently called the Goldberg equation

n = 20(i 2 + ij + j 2 ). (1)
The integers i and j in the Goldberg equation are in fact the components of a twodimensional Goldberg vector G = (i, j), sometimes also called Coxeter coordinates.

To avoid the mirror effect, we always assume that 0 ≤ i ≤ j and 0 < j. This vector determines the distance and positions of the vertices of the (i, j)triangle in the hexagonal lattice. See Figure 1 for a construction method of an (2, 3)triangle. Precisely 20 such (i, j)-triangles produce an (i, j)-icosahedral fullerene in a manner shown on Figure 2. The vertices of the triangles are centers of the 12 pentagons of the fullerene graph. The pair of the triangles ABC and EB ′ C ′ , as shown on Figure 2, is a pair of opposite triangles. There are exactly 10 such pairs. The pentagons with centers in A and E (Figure 2) are called antipodal pentagons (similarly, the pentagons with centers in B and B ′ , as well as those with centers C and C ′ are antipodal).

A C B G = (2, 3)
The icosahedral group I is the group of rotational symmetries of the icosahedron and dodecahedron, of order 60. This group is isomorphic to A 5 , the alternating group of even permutations of five objects. The full icosahedral group I h , is the point group of all symmetries of the icosahedron and dodecahedron. This group is equivalent to the group direct product A 5 × Z 2 of the alternating group A 5 and cyclic group Z 2 .

The order of the full icosahedral group of symmetries is 120. The (i, i)-and (0, i)- icosahedral fullerene graphs, i > 0, have full icosahedral symmetry group, i.e., every element of this class of graphs has a symmetry group I h .

A A A A A C B B D D E E E E E B ′ C ′
The idea of Goldberg was then used by Coxeter [START_REF] Coxeter | Virus macromolecules and geodesic domes[END_REF] and Graver [START_REF] Graver | Catalog of All Fullerenes with Ten or More Symmetries[END_REF] to classify highly symmetric fullerene graphs.

Nanotubical fullerene graphs

While the icosahedral fullerens have " spherical" shape, there is a class of fullerene graphs of tubular shapes, called nanotubical graphs or simply nanotubes. They are cylindrical in shape, with the two ends capped with a subgraph containing six pentagons and possibly some hexagons. The cylindrical part of the nanotube can be obtained by rolling a planar hexagonal grid. The way the grid is wrapped is represented by a Goldberg vector (i, j), called also the type of the nanotube. See Figure 3 for an example of the construction of the cylindrical part of a nanotube.

Nanotube fullerene graphs exist for all vectors (i, j) with i + j ≥ 6 and for the vector (0, 5). A (0, 5) nanotube is depicted on Figure 5.

Nanotubes with i = 0 are called zig-zag nanotubes, and the ones with i = j are called armchair nanotubes. These are the only types of nanotubes where the cylindrical part has a mirror symmetry. The Buckminsterfullerene C 60 can be viewed as the smallest nanotube of type [START_REF] Boltalina | Direct fluorination of fullerenes[END_REF][START_REF] Boltalina | Direct fluorination of fullerenes[END_REF], see Figure 4. It is the smallest nanotube with the caps satisfying the IPR.

There are such nanotubes for all types (i, j) with i + j ≥ 11 and for (5, 5), (0, 9), (0, 10), [START_REF] Andova | On the diameter and some related invariants of fullerene graphs[END_REF][START_REF] Došlić | On lower bounds of number of perfect matchings in fullerene graphs[END_REF], and (2, 8). = + + Figure 4: Buckminsterfullerene is the smallest nanotube of type [START_REF] Boltalina | Direct fluorination of fullerenes[END_REF][START_REF] Boltalina | Direct fluorination of fullerenes[END_REF].

Structural properties of fullerene graphs

As we already mentioned earlier, each fullerene graph has 12 pentagonal faces, all the other faces are hexagons. The distribution of the pentagons determines the shape of the corresponding polyhedron/molecule. If the pentagons are distributed uniformly, the polyhedron is spherical, with the pentagons being the source of the curvature.

Cyclic edge-cuts in fullerene graphs

Recall that an edge-cut

C of G is cyclic if each component of G -C contains a cycle.
A graph is cyclically k-edge-connected if at least k edges must be removed to disconnect it into two components such that each contains a cycle. For fullerene graphs, Došlić [START_REF] Došlić | Cyclical edge-connectivity of fullerene graphs and (k, 6)-cages[END_REF] proved the following theorem.

Theorem 2.1. Every fullerene graph F is cyclically 5-edge-connected.

Clearly, the result is best possible, since each pentagonal face is separated from the rest of the graph by 5 edges. We say that a cyclic 5-edge cut is trivial, if it separates only one pentagonal face. Kardoš and Škrekovski [START_REF] Kardoš | Cyclic edge-cuts in fullerene graphs[END_REF] proved that fullerenes which are not cyclically 6-edge-connected are of unique type.

Theorem 2.2. A (0, 5)-nanotube is the only fullerene graph with a nontrivial cyclic 5-edge-cut.

A (0, 5)-nanotube is depicted on Figure 5. Observe that every nontrivial cyclic 5-edge-cut of a (0, 5)-nanotube separates the two caps (containing 6 pentagonal faces each) from each other. Similarly, for each nanotube of type (i, j) there are cyclic (i+j)-edge cuts separating the two caps from each other. However, the nanotubes are not the only fullerene graphs with such edge-cuts. An example of a fullerene graph, in which the 12 pentagonal faces are partitioned into four clusters containing 1, 2, 4, and 5 pentagons each, is depicted on Figure 6. It is easy to see that this graph is not a nanotube.

Figure 6: A half of a fullerene graph, which is not a nanotube, but it is possible to separate six pentagonal face from the other six. To obtain the whole graph, the depicted graph is to be glued with its mirror copy along the dashed line.

Possible stability predictors

Fullerenes were studied from different perspective. The study from graph theoretical point of view has been motivated by a search for invariants that will correlate with their stability as a compound. Banhart et al. [START_REF] Banhart | The formation, annealing and self-comression of carbon onions under electron irradiation[END_REF] observed that molecular structures with sharp edges and facets are less stable than the structures with smooth curvature.

In this sense, the more molecule graph is spherical, the more it should be stable.

A number of graph-theoretical invariants were examined as potential stability predictors with various degrees of success [START_REF] Fajtlowicz | Graph-Theoretic Independence as a Predictor of Fullerene Stability[END_REF][START_REF] Došlić | Bipartivity of fullerene graphs and fullerene stability[END_REF]. Among the promising candidates for stability predictors are: the bipartite edge frustration, the independence number, the saturation number, the number of perfect matchings, etc. As a result of those investigations, we have achieved a fairly thorough understanding of fullerene graphs and their properties. However, some problems and questions still remain open [START_REF] Daugherty | Independent Sets and Closed-Shell Independent Sets of Fullerenes[END_REF][START_REF] Malkevitch | Geometrical and combinatorial questions about fullerenes[END_REF][START_REF] Fowler | Facts and conjectures about fullerene graphs: leapfrog, cylindrical and Ramanujan fullerenes[END_REF]. Special place among them have several interesting conjectures made by Graffiti, a conjecture making software [START_REF] Fajtlowicz | Written on the Wall (A List of Conjectures of Graffiti)[END_REF]. Here we consider several of those open questions starting from our results on the fullerene diameters and a recent result on the bipartite edge frustration [START_REF] Dvořak | Bipartizing fullerenes[END_REF][START_REF] Faria | Odd cycle transversals and independent sets in fullerene graphs[END_REF].

Diameter of fullerene graphs

For large enough fullerene graphs of spherical shape, the diameter diam(G) (the maximum length of the shortest path between two vertices) is proportional to the radius of the sphere, whereas the number of vertices n is proportional to its surface.

Hence, one could expect the diameter of such graphs to be of order Θ( √ n). On the other hand, for a nanotubical fullerene graph of type (i, j), the diameter diam(G) is proportional to the length of the cylindrical part of the graph, whereas the number of vertices n is proportional to the product of the length of the tube and its circumference i + j. In this case, one could expect the diameter of such graphs to be of order Θ(n).

A well known result on the degree-diameter problem states that the number of vertices in a planar graph with maximum degree 3 grows at most exponentially with diameter [START_REF] Friedman | New bounds for largest planar graphs with fixed maximum degree and diameter[END_REF].

Proposition 3.1. Let G be a planar graph with maximum degree 3. Then, G has at

most 2 diam(G)+1 -1 vertices.
This results gives a logarithmic lower bound on the diameter in terms of the number of vertices.

The logarithmic character of the bound can be attributed to the presence of faces of large size. It would be reasonable to expect that better lower bounds exist for polyhedral graphs with bounded face size.

As mentioned above, fullerene graphs only have pentagonal and hexagonal faces, and this fact can be used to show that the diameter is of order Ω( √ n) indeed.

The following theorem gives the lower bound for the diameter of fullerenes [START_REF] Andova | On the diameter and some related invariants of fullerene graphs[END_REF].

Theorem 3.2. Let G be a fullerene graph. Then,

diam(G) ≥ √ 24n -15 -3 6 .
For fullerene graphs with full icosahedral symmetry the diameter has been determined exactly in [START_REF] Andova | Diameter of full icosahedral-symmetry fullerene graphs[END_REF]:

Theorem 3.3. Let G be a (0, i)-icosahedral fullerene graph, with i > 0. Then diam(G) = 6i -1 = 9 5 n -1.
Theorem 3.4. Let G be an (i, i)-icosahedral fullerene graph, with i > 0. Then

diam(G) = 10i -1 = 5 3 n -1.
As we believe that the icosahedral fullerenes have the smallest diameter, it leads us to conjecture the following. As far as upper bound is concerned, the nanotubes of small circumference are the extremal graphs, as says the following theorem [START_REF] Andova | On the diameter and some related invariants of fullerene graphs[END_REF].

Theorem 3.5. Let F be a fullerene graph with n vertices. Then,

diam(F ) ≤ n 6 + 5 2 ,
unless F is a (0, 5)-nanotube. In that case, we have n = 10 k, k ≥ 3 and

diam(F ) = n 5 , k ∈ {3, 4} ; n 5 -1, k ≥ 5.

Bipartite edge frustration

The bipartite edge frustration of a graph G, denoted by ϕ(G), is the smallest cardinality of a set of edges of G that need to be removed from G in order to obtain a bipartite spanning subgraph. Bipartite edge frustration of a fullerene graph G can be efficiently computed by finding a minimum-weight perfect matching in the pentagondistance graph of G [START_REF] Došlić | Computing the bipartite edge frustration of fullerene graphs[END_REF]. In the same reference it was shown that ϕ(G) ≥ 6 for any fullerene graph G and that this bound is sharp. Furthermore, it was shown that the bipartite edge frustration of fullerene graphs with icosahedral symmetry is proportional to the square root of the number of vertices ( [START_REF] Došlić | Computing the bipartite edge frustration of fullerene graphs[END_REF], Proposition 11 and Corollary 12). The numerical computations suggested that it cannot behave worse than that, and prompted the authors to state the following conjecture.

Conjecture 3.2. Let G be a fullerene graph with n vertices. Then, ϕ(G) ≤ 12 5 n.

First, Dvořák, Lidický, and Škrekovski [START_REF] Dvořak | Bipartizing fullerenes[END_REF] proved a theorem with a weaker multiplicative constant. Later Faria, Klein and Stehlík proved the following theorem, which implies Conjecture 3.2 [START_REF] Faria | Odd cycle transversals and independent sets in fullerene graphs[END_REF].

Theorem 3.6. Let G be a fullerene graph with n vertices. Then, ϕ(G) ≤ 12 5 n. Moreover, the equality holds if and only if G is an (i, i)-icosahedral fullerene graph.

The bipartite edge frustration is a very reliable predictor of molecular stability: If ϕ(G) is close to the upper bound, it means that there is no way to partition the 12 pentagonal faces into six pairs in such a way that the pentagons within the pairs are close to each other.

On the other way, a small value for ϕ(G) does not say much about the structure of the fullerene graph and about the shape of the molecule. The smallest possible value of ϕ(G) is attained for any fullerene graph where the 12 pentagonal faces can be partitioned into 6 pairs, adjacent to each other. Those pairs of pentagonal faces can be distributed in a widely varied list of possible configurations, leading to molecules of different shapes, including nanotubes of all possible types.

Independence number

Another invariant tested as a possible stability predictor is the independence number [START_REF] Fajtlowicz | Graph-Theoretic Independence as a Predictor of Fullerene Stability[END_REF]. A set I ⊆ V (G) is independent if no two vertices from I are adjacent in G. The cardinality of any largest independent set in G is called the independence number of G and denoted by α(G). Independence number of fullerene graphs attracted a lot of attention also in the context of study of independent sets as possible models for addition of bulky segregated groups such as free radicals or halogen atoms [START_REF] Boltalina | Direct fluorination of fullerenes[END_REF].

Sharp upper bounds on the independence number of n/2 -2 for general fullerenes and n/2 -4 for those with isolated pentagons follow by simple counting argument [START_REF] Fowler | Independence number and fullerene stability[END_REF]. Lower bounds were gradually improved from (almost) trivial α(G) ≥ n/3 valid for all 3-chromatic graphs to α(G) ≥ 3 8 n [START_REF] Heckman | Independent sets in triangle-free cubic planar graphs[END_REF] valid for all triangle-free planar cubic graphs. A better lower bound of type α(G) ≥ n 2 -C √ n, for some constant C, was first established for icosahedral fullerenes [START_REF] Graver | The independence number of fullerenes and benzenoids[END_REF]. The first result that approaches the conjecture is due to Andova et al. [START_REF] Andova | On the diameter and some related invariants of fullerene graphs[END_REF].

Theorem 3.7. Let G be a fullerene graph with n vertices. Then

α(G) ≥ n 2 -78.58 √ n .
Finally, using Theorem 3.6, Conjecture 3.3 was confirmed it the following form.

Theorem 3.8. Let G be a fullerene graph with n vertices. Then,

α(G) ≥ n 2 - 3 5 n .
Moreover, the equality holds if and only if G is an (i, i)-icosahedral fullerene graph.

Similarly to the bipartite edge frustration, a small value of independence number means that the corresponding polyhedron/molecule is close to being of a spherical shape. On the other hand, a value of α(G) close to n/2 does not say much about the structure or the shape of the molecule.

The relations between diameter and the independence number of fullerenes appear in Conjecture 912 of Graffiti [START_REF] Fowler | Facts and conjectures about fullerene graphs: leapfrog, cylindrical and Ramanujan fullerenes[END_REF]. This conjecture was established in [START_REF] Andova | On the diameter and some related invariants of fullerene graphs[END_REF] for large fullerenes using Theorems 3.5 and 3.7.

Conjecture 3.4. If G is a fullerene graph, then α(G) ≥ 2(diam(G) -1).
Finally, Theorem 3.5 and Theorem 3.8 together confirm the conjecture in the affirmative.

Upper bound on the smallest eigenvalue

Another invariant investigated for its stability-predicting potential is the smallest eigenvalue of a fullerene graph. In [START_REF] Fowler | A note on the smallest eigenvalue of fullerenes[END_REF] it was proven that the dodecahedron has maximum smallest eigenvalue among all fullerene graphs and it is equal to -√ 5. It was also shown that the buckminsterfullerene C 60 has the maximum smallest eigenvalue among all IPR fullerene graphs. This observation lead the authors to state the following conjecture on the smallest eigenvalue of fullerene graphs with at least 60 vertices.

Conjecture 3.5. Among all fullerene graphs with at least 60 vertices, the buckminster fullerene has the maximum smallest eigenvalue.

The Conjectures 3.

λ n (G) = 3 -µ ∞ (G).
By plugging this into Theorem 3.9 and noting that bip(G) =

Number of perfect matchings

Since all carbon atoms are 4-valent, for every atom precisely one of the three bonds should be doubled. Such a set of double bonds is called a Kekulé structure in a fullerene. It corresponds to the notion of perfect matchings in fullerene graphs: a matching in a graph G is a set of edges of G such that no two edges in M share an end-vertex. A matching M is perfect if any vertex of G is incident with an edge of M.

Let M be a perfect matching in a fullerene graph G. A hexagonal face is resonant if it is incident with three edges in M. The maximum size of a set of resonant hexagons in G is called the Clar number of G.

The computation of the number of perfect matchings in typical fullerene graphs with a small number of vertices [START_REF] Došlić | Fullerene graphs with exponentially many perfect matchings[END_REF] indicates that this number should grow exponentially with n.

The first general lower bounds for the number of perfect matchings in fullerene graphs were linear in the number of vertices [START_REF] Došlić | On lower bounds of number of perfect matchings in fullerene graphs[END_REF][START_REF] Došlić | On some structural properties of fullerene graphs[END_REF][START_REF] Zhang | New lower bound on the number of perfect matchings in fullerene graphs[END_REF]. The best of them asserts that a fullerene graph with n vertices has at least 3(n+2) 4 different perfect matchings [START_REF] Zhang | New lower bound on the number of perfect matchings in fullerene graphs[END_REF].

The number of perfect matchings was proved to be exponential for several special classes of fullerene graphs: (0, 5)-nanotubes [START_REF] Kutnar | On cyclic edge-connectivity of fullerenes[END_REF], icosahedral fullerenes [START_REF] Došlić | Fullerene graphs with exponentially many perfect matchings[END_REF] or those who are obtained using specific operations [START_REF] Došlić | Leapfrog fullerenes have many perfect matchings[END_REF], before the following general result [START_REF] Kardoš | Fullerene graphs have exponentially many perfect matchings[END_REF]. The idea of the proof is to find a perfect matching with (n -380)/61 non-adjacent resonant hexagons. For such a perfect matching, any subset of the set of resonant hexagons can be switched to obtain a different perfect matching.

The bound from Theorem 3.11 is far from being optimal. A little is known about the relationship between the number of perfect matchings of a fullerene graph and the shape of the corresponding molecule.

Saturation number

The last invariant considered here, the saturation number, is also related to matchings.

The existence of perfect (and hence maximum) matchings in fullerene graphs has been established long time ago, and there are many papers concerned with their structural and enumerative properties [START_REF] Došlić | On some structural properties of fullerene graphs[END_REF][START_REF] Kardoš | Fullerene graphs have exponentially many perfect matchings[END_REF][START_REF] Klein | Theorems on carbon cages[END_REF]. Another class of matchings, the maximal matchings, have received much less attention so far, in spite of being potentially useful as mathematical models of dimer absorption. A matching M is maximal if it cannot be extended to a larger matching of G. The saturation number of G, denoted by s(G), is the cardinality of any smallest maximal matching of G. The saturation number of fullerene graph was studied in [START_REF] Došlić | On some structural properties of fullerene graphs[END_REF][START_REF] Došlić | Saturation number of fullerene graphs[END_REF], where the following bounds were established.

Theorem 3.12. There exists an absolute constant C such that

3n 10 ≤ s(G) ≤ n 2 -C log 2 n ,
for any fullerene graph G with n vertices.

Using the lower bound on the diameter we are able to improve the logarithmic additive correction of Theorem 3.12 and to prove the following result [START_REF] Andova | On the diameter and some related invariants of fullerene graphs[END_REF]. Still, this is not the best possible bound of the saturation number of the fullerene graphs. In [START_REF] Andova | Bounds on saturation number on fullerene graphs[END_REF] we prove that the saturation number for fullerenes on n vertices is essentially n/3.

Theorem 3.14. Let F be a fullerene graph on n vertices. Then,

n 3 -2 ≤ s(F ) ≤ n 3 + O( √ n) .
In order to prove the lower bound of this theorem we used the discharging method.

For the upper bound we first used the Theorem 3.6, and obtained a bipartite graph F ′ . Later, we establish that F ′ is an induced subgraph of a hexagonal lattice, or an induced subgraph of a hexagonal tube (defined as on Figure 3). Then define a maximal matching on F ′ such that from each hexagon precisely four vertices are covered by the matching.

Figure 1 :

 1 Figure 1: Construction of a (2, 3)-triangle, a metaface of a (2, 3)-icosahedral fullerene graph. The vertices of the equilateral triangle ABC are centers of pentagons.

Figure 2 :

 2 Figure 2: A (2, 3)-icosahedral fullerene. Its triangular faces are constructed as on Figure 1. The vertices with a same name coincide. The vertices of each triangular are centers of the pentagons. The pentagons with center A and center E are antipodal pentagons, and the triangles ABC and EB ′ C ′ form a pair of opposite triangles.

Figure 3 :

 3 Figure 3: An example of a (2, 4) nanotube. The hexagons denoted equally overlap.

Figure 5 :

 5 Figure 5: A (0, 5)-nanotube.

Conjecture 3 . 1 .

 31 Let G be a fullerene graph on n vertices. Then diam(

Conjecture 3 . 3 . 2 -

 332 This observation was formalized in a pair of conjectures in a recent Ph.D. thesis by S.Daugherty ([8, pp. 96]). The first one states that the minimum possible independence number is achieved on the icosahedral fullerenes that also figure prominently in Conjecture 3.2. The second one [8, Conjecture 5.5.2] states the precise form of the conjectured lower bound. Notice that the constant 3/ √ 15 is exactly one half of the constant 12/5 of Conjecture 3.2. Let G be a fullerene graph with n vertices. Then α(G) ≥ n 3 n/15.

Theorem 3 . 9 .

 39 2 and 3.5 are connected via a result on Laplacian eigenvalues from the monograph by Godsil and Royle ([28, pp. 293]). Let G be a graph with n vertices. Then bip(G) ≤ n 4 µ ∞ (G). Here bip(G) denotes the maximum number of edges in a bipartite spanning subgraph of G (hence the number of edges in G minus the bipartite edge frustration), and µ ∞ (G) is the largest Laplacian eigenvalue of G. Recall that the smallest eigenvalue λ n (G) of a 3-regular graph G and the largest Laplacian eigenvalue µ ∞ (G) of G are related via the following relation ([28, pp. 280]):

3 2 n

 2 -ϕ(G) we obtain an upper bound on λ n (G) in terms of the bipartite edge frustration of G of the form λ n (G) ≤ -3 + 4 n ϕ(G). By taking into account an upper bound on ϕ(G) we immediately obtain the following upper bound on the smallest eigenvalue of a fullerene graph with n vertices [20]. Theorem 3.10. Let G be a fullerene graph with n vertices. Then, λ n (G) ≤ -3

Theorem 3 . 11 .

 311 Let G be a fullerene graph with n vertices that has no non-trivial cyclic 5-edge-cut. The number of perfect matchings of G is at least 2 n-380 61 .

Theorem 3 . 13 .

 313 Let G be a fullerene graph with n vertices. Then,

Partially supported by Slovenian ARRS Program P1-00383 and Creative Core -FISNM -3330-13-500033.