- Email: [email protected]
Rigorous representation and exact simulation of real Gaussian stationary processes
Chemical Physics 375 (2010) 378–379
Contents lists available at ScienceDirect
Chemical Physics journal homepage: www.elsevier.com/locate/chemphys
Rigorous representation and exact simulation of real Gaussian stationary processes Jiushu Shao * College of Chemistry, Beijing Normal University, Beijing 100875, China
a r t i c l e
i n f o
Article history: Available online 25 June 2010 Keywords: Gaussian process Numerical simulation
a b s t r a c t It is shown that the real Gaussian stationary process with zero mean can be exactly interpreted as a convolution defined by its covariance and a complex white noise plus the conjugate pair of the complex white noise. When the algorithm for generating the white noise is given, such a representation thereby offers a direct and convenient technique to numerically simulate Gaussian processes with arbitrary correlations. Moreover, the fast Fourier transform can be invoked for solving the convolution, which leads to a feasible implementation of the simulation method. As a primary test, a Gaussian process in quantum dissipation is simulated. This stochastic process has an infinitely long correlation time and is hardly treated by other approaches. Ó 2010 Elsevier B.V. All rights reserved.
Due to the unavoidably intrinsic or external uncertainty, the evolution of a macroscopic system often manifests as stochastic processes [1–3]. The best known example is the Brownian motion, which characterizes the dynamical behavior of a particle subjected to random forces resulting from the environment or heat bath of many degrees of freedom [4]. It is observed that the Brownian motion can be described by a stochastic differential equation, the socalled Langevin equation. In principle, one may derive the equation of motion for a dissipative system from a microscopic model and thus be able to solve the dynamics once the initial condition and the random force posed by the heat bath are specified [5]. The quantum counterpart of dissipative dynamics shares the same spirit, namely, the density matrix of the system evolves in the stochastic field induced by the heat bath [6]. Compared to the noisefree case, the stochastic dynamics often displays diversified features and the role of the stochastic field in the dynamics is also versatile. Of most interesting is that noises can be counterintuitively constructive. A good example is the stochastic resonance [7], an effect that might have been used in the biological world [8]. To understand this and other noise-induced effects one has to solve and analyze the underlying equation of motion. There are many studies on the stochastic differential equation and how to solve it, from both theoretical and practical points of view [9,10]. As one can imagine, an indispensable step in simulating stochastic dynamics is to generate the involved random field and the computational cost of the simulation heavily depends on the feasibility of the generation of the random field. This paper addresses the issue of generating stochastic fields with given statistical properties. To be specific, we propose a rep-
* Tel.: +86 10 58802417; fax: +86 10 58802075. E-mail address: [email protected] 0301-0104/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.chemphys.2010.06.027
resentation of a real Gaussian stationary process in terms of its covariance and complex white noises. This formulation allows us to develop a computationally feasible algorithm for generating any Gaussian stochastic processes with given averages and covariances. Traditionally, there are several ways to simulate Gaussian stochastic process [11], which are exemplified by the Cholesky decomposition and Durbin–Levinson (DL) recursion [12,13], the Davies–Harte (DH) or circulant embedding approach [14–16], or methods based on the spectral analysis [17]. The DL algorithm is valid for any reasonable covariance, but its computational cost scales as O(N2) for a process of N time steps. By contrast, the DH scheme scales as O(Nlog(N)) in numerical implementation. Although the DH method is of superior performance and thus gains many applications [15,16,18–20], its validity is severely limited by the requirement of nonnegative covariance or negative covariance for nonzero time lags. Let us think of the numerical simulation for a real Gaussian stationary process v(t) defined by the mean hv(t)i that is assumed to be zero without loss of generality and the covariance hv(t)v(t0 )i = a(tt0 ). Now consider a complex stochastic process described by
nðtÞ ¼
Z 0
t
1 0 dt aðt t 0 Þ½l1 ðt0 Þ þ il2 ðt 0 Þ þ ½l1 ðtÞ il2 ðtÞ; 2
ð1Þ
where l1,2(t) are two real, independent white noises, hl1,2(t)i = 0, hlj(t)lj(t0 )i = d(t t0 ) (j = 1, 2). It is straightforward to show that both n(t) and its complex conjugate n*(t) satisfy the same statistical properties as v(t) and thus the former can be regarded as a representation of the latter. The equivalence of n(t) and v(t) was recognized in the stochastic description of quantum dissipative systems [21] although the relationship was not then used as a tool for generating Gaussian processes. This formulation is not limited to one dimension. In fact, for a real multidimensional stochastic process
J. Shao / Chemical Physics 375 (2010) 378–379
379
~ wðtÞ with covariance matrix C(t) one can establish a similar representation ~ /ðtÞ, namely
Z
~ /ðtÞ ¼
t
0
1 0 dt Cðt t0 Þ½~ l1 ðt0 Þ þ i~ l2 ðt0 Þ þ ½~ l1 ðtÞ i~ l2 ðtÞ; 2
ð2Þ
where ~ l1;2 ðtÞ are multidimensional white noises. For brevity, only one-dimensional case will be considered in the following. Recognizing that the convolution integral in Eq. (1) can be calculated via fast Fourier transform (FFT), one immediately works out an efficient means of simulating the real Gaussian process with arbitrary covariance. In particular, we want to create a time series n(t) at 0, Dt, . . . , NDt with Dt = t/N, which is a realization of the stochastic process v(t). Allowing for a good simulation even for the zero-time covariance, one may recast Eq. (1) in the following discrete form:
nj ¼ Dt
j1 X
1 2
1 2
ajk ðl1k þ il2k Þ þ Dta0 ðl1j þ il2j Þ þ ðl1j il2j Þ:
k¼0
ð3Þ Generating nj for j = 0 to j = N thus requires a white noise generator and fast discrete convolution manipulation through the FFT algorithm. To this end we can use several well known programs in literature say from ”Numerical Recipes” [22]. Note that for calculating the convolution one should employ 2N-point instead of Npoint FFT in which aj, l1j and l2j are zero-padded for j > N. There remains flexibility in the representation of the real stochastic process v(t), Eq. (1). For instance, the upper bound in the convolution can be changed to any value tu as long as tu > t. Moreover, Eq. (1) can be modified upon introducing a finite, yet arbitrary function s(t) such that
nðtÞ ¼ sðtÞ
Z 0
t
0
dt aðt t0 Þ½l1 ðt 0 Þ þ il2 ðt 0 Þ þ
1 ½l ðtÞ il2 ðtÞ: 2sðtÞ 1 ð4Þ
Therefore, it is free to choose a scaling function s(t) or simply a constant to make the simulation of the stochastic process v (t) more effective. By effective we mean that fewer samplings are needed to reach the stochastic average of v(t). Let us now implement the approach to simulate the stochastic field in the spin-boson model [23,24]. At zero temperature, the Ohmic dissipation corresponds exactly to two kinds of stochastic Gaussian fields, one with a very short memory and the other a long one. The short-memory noise is readily treated by virtue of deterministic method [21,24] and here the long-range one is interesting to us. With an exponential cutoff its covariance assumes the form
aðtÞ ¼
1 t2 ð1 þ t 2 Þ2
:
When t is large, the stochastic process defies a fast simulation via DH method. But this is not a problem in our approach. Fig. 1. displays the results of the calculations and clearly shows that this scheme works. The covariance of the generated time series is estiPM1 ðkÞ ðkÞ mated through aj ¼ k¼0 n1 nj =M where k dictates the kth realization and M the total number of samplings. In the simulations the scaling function s(t) is taken to be a constant 9. To summarize, we propose a rigorous representation for a real Gaussian stationary process. Based on this description, a feasible simulation algorithm is developed. Compared to other fast generators of stochastic processes, no restriction on the covariance is required. This approach is tested through an example resulting from quantum dissipation and promising results are obtained. Nevertheless, there is a drawback of this method, namely, the replacement of the real stochastic process replaced by a complex substitute, which may cause ineffectiveness in some applications.
Fig. 1. Covariance: Comparison of the exact one with the simulation results. (a) is from the exact calculation. (b) is from the estimation of 4 104 realizations. (c) is from the estimation of 16 104 realizations.
Acknowledgment The paper is dedicated to Professor Peter Hänggi, a great scientist and mentor on the occasion of his 60th birthday. This work is supported by the National Natural Science Foundation of China (Nos. 20773017 and 20533060) and the 973 program of the Ministry of Science and Technology of China (2007CB815206). References [1] N.G. van Kampen, Stochastic Processes in Physics and Chemistry, second ed., Elsevier, Amsterdam, 1997. [2] P. Hänggi, H. Thomas, Phys. Rep. 88 (1982) 207. [3] C.W. Gardiner, Handbook of Stochastic Methods, third ed., Springer, Berlin, 2004. [4] R. Kubo, M. Toda, N. Hashitume, Statistical Physics II, second ed., Springer, Berlin, 1992. [5] R. Zwanzig, Nonequilibrium Statistical Mechanics, Oxford Univ Press, Oxford, 2001. [6] R. Kubo, J. Math. Phys. 4 (1963) 174; J.T. Stockburger, H. Grabert, Phys. Rev. Lett. 88 (2002) 170407; J. Shao, J. Chem. Phys. 120 (2004) 5053. [7] L. Gammaitoni, P. Hänggi, P. Jung, F. Marchesoni, Rev. Mod. Phys. 70 (1998) 223; T. Wellens, V. Shatokhin, A. Buchleitner, Rep. Prog. Phys. 67 (2004) 45. [8] P. Hänggi, Chem. Phys. Chem. 3 (2002) 285; M.D. McDonnell, D. Abbott, PLoS Comp. Bio. 5 (2009) e1000348. [9] B. Øksendal, Stochastic Differential Equations, sixth ed., Springer, Berlin, 2005; P.E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer, Berlin, 2000. [10] D.T. Gillespie, J. Phys. Chem. 81 (1977) 2340; D.T. Gillespie, J. Comput. Phys. 22 (1976) 403; D.T. Gillespie, J. Chem. Phys. 113 (2000) 297; C. Kim, E.K. Lee, P. Hänggi, Peter Talkner, Phys. Rev. E 76 (2007) 011109. [11] G.E. Johnson, Proc. IEEE 82 (1994) 270. [12] C.E. Rasmussen, C.K.I. Williams, Gaussian Processes for Machine Learning, MIT Press, Cambridge, 2006. [13] P.J. Brockwell, R.A. Davis, Time Series: Theory and Methods, second ed., Springer, New York, 1991. [14] R.B. Davies, D.S. Harte, Biometrika 74 (1987) 95. [15] C.R. Dietrich, G.N. Newsam, Wat. Resour. Res. 29 (1993) 2861. [16] A.T.A. Wood, G. Chan, J. Comput. Graph. Stat. 3 (1994) 409. [17] H.A. Makse, S. Havlin, M. Schwartz, H.E. Stanley, Phys. Rev. E 53 (1996) 5445. [18] P.F. Craigmile, J. Time Series Anal. 24 (2003) 505. [19] D.B. Percival, Signal Proc. 86 (2006) 1470. [20] T. Gneiting, H. Ševcˇíková, D.B. Percival, M. Schlather, Y. Jiang, J. Comput. Graph. Stat. 15 (2006) 483. [21] Y. Zhou, Y. Yan, J. Shao, Europhys. Lett. 72 (2005) 334. [22] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes, third ed., Cambridge Univ Press, 2007. [23] A.J. Leggett, S. Chakravarty, A.T. Dorsey, M.P.A. Fisher, A. Garg, W. Zwerger, Rev. Mod. Phys. 59 (1987) 1. [24] Y. Zhou, J. Shao, J. Chem. Phys. 128 (2008) 034106.
Contents lists available at ScienceDirect
Chemical Physics journal homepage: www.elsevier.com/locate/chemphys
Rigorous representation and exact simulation of real Gaussian stationary processes Jiushu Shao * College of Chemistry, Beijing Normal University, Beijing 100875, China
a r t i c l e
i n f o
Article history: Available online 25 June 2010 Keywords: Gaussian process Numerical simulation
a b s t r a c t It is shown that the real Gaussian stationary process with zero mean can be exactly interpreted as a convolution defined by its covariance and a complex white noise plus the conjugate pair of the complex white noise. When the algorithm for generating the white noise is given, such a representation thereby offers a direct and convenient technique to numerically simulate Gaussian processes with arbitrary correlations. Moreover, the fast Fourier transform can be invoked for solving the convolution, which leads to a feasible implementation of the simulation method. As a primary test, a Gaussian process in quantum dissipation is simulated. This stochastic process has an infinitely long correlation time and is hardly treated by other approaches. Ó 2010 Elsevier B.V. All rights reserved.
Due to the unavoidably intrinsic or external uncertainty, the evolution of a macroscopic system often manifests as stochastic processes [1–3]. The best known example is the Brownian motion, which characterizes the dynamical behavior of a particle subjected to random forces resulting from the environment or heat bath of many degrees of freedom [4]. It is observed that the Brownian motion can be described by a stochastic differential equation, the socalled Langevin equation. In principle, one may derive the equation of motion for a dissipative system from a microscopic model and thus be able to solve the dynamics once the initial condition and the random force posed by the heat bath are specified [5]. The quantum counterpart of dissipative dynamics shares the same spirit, namely, the density matrix of the system evolves in the stochastic field induced by the heat bath [6]. Compared to the noisefree case, the stochastic dynamics often displays diversified features and the role of the stochastic field in the dynamics is also versatile. Of most interesting is that noises can be counterintuitively constructive. A good example is the stochastic resonance [7], an effect that might have been used in the biological world [8]. To understand this and other noise-induced effects one has to solve and analyze the underlying equation of motion. There are many studies on the stochastic differential equation and how to solve it, from both theoretical and practical points of view [9,10]. As one can imagine, an indispensable step in simulating stochastic dynamics is to generate the involved random field and the computational cost of the simulation heavily depends on the feasibility of the generation of the random field. This paper addresses the issue of generating stochastic fields with given statistical properties. To be specific, we propose a rep-
* Tel.: +86 10 58802417; fax: +86 10 58802075. E-mail address: [email protected] 0301-0104/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.chemphys.2010.06.027
resentation of a real Gaussian stationary process in terms of its covariance and complex white noises. This formulation allows us to develop a computationally feasible algorithm for generating any Gaussian stochastic processes with given averages and covariances. Traditionally, there are several ways to simulate Gaussian stochastic process [11], which are exemplified by the Cholesky decomposition and Durbin–Levinson (DL) recursion [12,13], the Davies–Harte (DH) or circulant embedding approach [14–16], or methods based on the spectral analysis [17]. The DL algorithm is valid for any reasonable covariance, but its computational cost scales as O(N2) for a process of N time steps. By contrast, the DH scheme scales as O(Nlog(N)) in numerical implementation. Although the DH method is of superior performance and thus gains many applications [15,16,18–20], its validity is severely limited by the requirement of nonnegative covariance or negative covariance for nonzero time lags. Let us think of the numerical simulation for a real Gaussian stationary process v(t) defined by the mean hv(t)i that is assumed to be zero without loss of generality and the covariance hv(t)v(t0 )i = a(tt0 ). Now consider a complex stochastic process described by
nðtÞ ¼
Z 0
t
1 0 dt aðt t 0 Þ½l1 ðt0 Þ þ il2 ðt 0 Þ þ ½l1 ðtÞ il2 ðtÞ; 2
ð1Þ
where l1,2(t) are two real, independent white noises, hl1,2(t)i = 0, hlj(t)lj(t0 )i = d(t t0 ) (j = 1, 2). It is straightforward to show that both n(t) and its complex conjugate n*(t) satisfy the same statistical properties as v(t) and thus the former can be regarded as a representation of the latter. The equivalence of n(t) and v(t) was recognized in the stochastic description of quantum dissipative systems [21] although the relationship was not then used as a tool for generating Gaussian processes. This formulation is not limited to one dimension. In fact, for a real multidimensional stochastic process
J. Shao / Chemical Physics 375 (2010) 378–379
379
~ wðtÞ with covariance matrix C(t) one can establish a similar representation ~ /ðtÞ, namely
Z
~ /ðtÞ ¼
t
0
1 0 dt Cðt t0 Þ½~ l1 ðt0 Þ þ i~ l2 ðt0 Þ þ ½~ l1 ðtÞ i~ l2 ðtÞ; 2
ð2Þ
where ~ l1;2 ðtÞ are multidimensional white noises. For brevity, only one-dimensional case will be considered in the following. Recognizing that the convolution integral in Eq. (1) can be calculated via fast Fourier transform (FFT), one immediately works out an efficient means of simulating the real Gaussian process with arbitrary covariance. In particular, we want to create a time series n(t) at 0, Dt, . . . , NDt with Dt = t/N, which is a realization of the stochastic process v(t). Allowing for a good simulation even for the zero-time covariance, one may recast Eq. (1) in the following discrete form:
nj ¼ Dt
j1 X
1 2
1 2
ajk ðl1k þ il2k Þ þ Dta0 ðl1j þ il2j Þ þ ðl1j il2j Þ:
k¼0
ð3Þ Generating nj for j = 0 to j = N thus requires a white noise generator and fast discrete convolution manipulation through the FFT algorithm. To this end we can use several well known programs in literature say from ”Numerical Recipes” [22]. Note that for calculating the convolution one should employ 2N-point instead of Npoint FFT in which aj, l1j and l2j are zero-padded for j > N. There remains flexibility in the representation of the real stochastic process v(t), Eq. (1). For instance, the upper bound in the convolution can be changed to any value tu as long as tu > t. Moreover, Eq. (1) can be modified upon introducing a finite, yet arbitrary function s(t) such that
nðtÞ ¼ sðtÞ
Z 0
t
0
dt aðt t0 Þ½l1 ðt 0 Þ þ il2 ðt 0 Þ þ
1 ½l ðtÞ il2 ðtÞ: 2sðtÞ 1 ð4Þ
Therefore, it is free to choose a scaling function s(t) or simply a constant to make the simulation of the stochastic process v (t) more effective. By effective we mean that fewer samplings are needed to reach the stochastic average of v(t). Let us now implement the approach to simulate the stochastic field in the spin-boson model [23,24]. At zero temperature, the Ohmic dissipation corresponds exactly to two kinds of stochastic Gaussian fields, one with a very short memory and the other a long one. The short-memory noise is readily treated by virtue of deterministic method [21,24] and here the long-range one is interesting to us. With an exponential cutoff its covariance assumes the form
aðtÞ ¼
1 t2 ð1 þ t 2 Þ2
:
When t is large, the stochastic process defies a fast simulation via DH method. But this is not a problem in our approach. Fig. 1. displays the results of the calculations and clearly shows that this scheme works. The covariance of the generated time series is estiPM1 ðkÞ ðkÞ mated through aj ¼ k¼0 n1 nj =M where k dictates the kth realization and M the total number of samplings. In the simulations the scaling function s(t) is taken to be a constant 9. To summarize, we propose a rigorous representation for a real Gaussian stationary process. Based on this description, a feasible simulation algorithm is developed. Compared to other fast generators of stochastic processes, no restriction on the covariance is required. This approach is tested through an example resulting from quantum dissipation and promising results are obtained. Nevertheless, there is a drawback of this method, namely, the replacement of the real stochastic process replaced by a complex substitute, which may cause ineffectiveness in some applications.
Fig. 1. Covariance: Comparison of the exact one with the simulation results. (a) is from the exact calculation. (b) is from the estimation of 4 104 realizations. (c) is from the estimation of 16 104 realizations.
Acknowledgment The paper is dedicated to Professor Peter Hänggi, a great scientist and mentor on the occasion of his 60th birthday. This work is supported by the National Natural Science Foundation of China (Nos. 20773017 and 20533060) and the 973 program of the Ministry of Science and Technology of China (2007CB815206). References [1] N.G. van Kampen, Stochastic Processes in Physics and Chemistry, second ed., Elsevier, Amsterdam, 1997. [2] P. Hänggi, H. Thomas, Phys. Rep. 88 (1982) 207. [3] C.W. Gardiner, Handbook of Stochastic Methods, third ed., Springer, Berlin, 2004. [4] R. Kubo, M. Toda, N. Hashitume, Statistical Physics II, second ed., Springer, Berlin, 1992. [5] R. Zwanzig, Nonequilibrium Statistical Mechanics, Oxford Univ Press, Oxford, 2001. [6] R. Kubo, J. Math. Phys. 4 (1963) 174; J.T. Stockburger, H. Grabert, Phys. Rev. Lett. 88 (2002) 170407; J. Shao, J. Chem. Phys. 120 (2004) 5053. [7] L. Gammaitoni, P. Hänggi, P. Jung, F. Marchesoni, Rev. Mod. Phys. 70 (1998) 223; T. Wellens, V. Shatokhin, A. Buchleitner, Rep. Prog. Phys. 67 (2004) 45. [8] P. Hänggi, Chem. Phys. Chem. 3 (2002) 285; M.D. McDonnell, D. Abbott, PLoS Comp. Bio. 5 (2009) e1000348. [9] B. Øksendal, Stochastic Differential Equations, sixth ed., Springer, Berlin, 2005; P.E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer, Berlin, 2000. [10] D.T. Gillespie, J. Phys. Chem. 81 (1977) 2340; D.T. Gillespie, J. Comput. Phys. 22 (1976) 403; D.T. Gillespie, J. Chem. Phys. 113 (2000) 297; C. Kim, E.K. Lee, P. Hänggi, Peter Talkner, Phys. Rev. E 76 (2007) 011109. [11] G.E. Johnson, Proc. IEEE 82 (1994) 270. [12] C.E. Rasmussen, C.K.I. Williams, Gaussian Processes for Machine Learning, MIT Press, Cambridge, 2006. [13] P.J. Brockwell, R.A. Davis, Time Series: Theory and Methods, second ed., Springer, New York, 1991. [14] R.B. Davies, D.S. Harte, Biometrika 74 (1987) 95. [15] C.R. Dietrich, G.N. Newsam, Wat. Resour. Res. 29 (1993) 2861. [16] A.T.A. Wood, G. Chan, J. Comput. Graph. Stat. 3 (1994) 409. [17] H.A. Makse, S. Havlin, M. Schwartz, H.E. Stanley, Phys. Rev. E 53 (1996) 5445. [18] P.F. Craigmile, J. Time Series Anal. 24 (2003) 505. [19] D.B. Percival, Signal Proc. 86 (2006) 1470. [20] T. Gneiting, H. Ševcˇíková, D.B. Percival, M. Schlather, Y. Jiang, J. Comput. Graph. Stat. 15 (2006) 483. [21] Y. Zhou, Y. Yan, J. Shao, Europhys. Lett. 72 (2005) 334. [22] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes, third ed., Cambridge Univ Press, 2007. [23] A.J. Leggett, S. Chakravarty, A.T. Dorsey, M.P.A. Fisher, A. Garg, W. Zwerger, Rev. Mod. Phys. 59 (1987) 1. [24] Y. Zhou, J. Shao, J. Chem. Phys. 128 (2008) 034106.