Eigenfunctions of the MultidimensionalLinear Noise Fokker-Planck Operator via Ladder Operators

# Eigenfunctions of the Multidimensional Linear Noise Fokker-Planck Operator via Ladder Operators

Todd K. Leen111Todd.Leen@georgetown.edu Graduate School of Arts and Sciences, Georgetown University Robert Friel Courant Institute of Mathematical Sciences, New York, NY David Nielsen
###### Abstract

The eigenfunctions and eigenvalues of the Fokker-Planck operator with linear drift and constant diffusion are required for expanding time-dependent solutions and for evaluating our recent perturbation expansion for probability densities governed by a nonlinear master equation. Although well-known in one dimension, for multiple dimensions the eigenfunctions are not explicitly given in the literature. We develop raising and lowering operators for the Fokker-Planck (FP) operator and its adjoint, and use them to obtain expressions for the corresponding eigenvalues and eigenfunctions. We show that the eigenfunctions for the forward and adjoint FP operators form a bi-orthogonal set, and that the eigenfunctions reduce to sums of products of Hermite functions in a particular coordinate system.

Keywords: Ornstein-Uhlenbeck operator, Fokker-Planck eigenfunctions, Hermite functions, ladder operators.

## 1 Introduction

The Fokker-Planck equation (FPE) with linear drift and constant diffusion describes an Ornstein-Uhlenbeck process. In one dimension, the eigenfunctions are the well-known Hermite functions (Risken, 1989; Gardiner, 2009). The eigenfunctions enable expansion of time-dependent solutions to the FPE, and are required to evaluate our recent perturbation expansion for the probability densities arising from a nonlinear master equation (Leen and Friel, 2011, 2012; Leen et al., 2012). (Thomas and Grima have recently derived a similar approximation and applied it to one-dimensional chemical and gene expression systems (2015).) The multi-dimensional eigenfunctions are not given in convenient form in the literature. (Liberzon and Brockett (2000) discuss the eigenvalue spectrum, but do not give an explicit form of the complete set of eigenfunctions, nor discuss the eigenfunctions of the adjunct operator.) In the restricted case that the drift and diffusion are simultaneously diagonalizable, the equations separate and the solution reduces to products of the one-dimensional eigenfunctions. However for general drift and diffusion, the eigenfunctions are not simply obtained from the differential operators.

This note gives raising and lowering operators for the forward and adjoint Fokker-Planck operators that develop the corresponding complete sets of eigenfunctions. We give bi-orthogonality relations and show that in a special coordinate system, the eigenfunctions are sums of products of Hermite functions.

## 2 Operators and Eigenfunctions

Let , , . The forward and adjoint (backward) Fokker-Planck operators for an N-dimensional Ornstein-Uhlenbeck process are

 Lq(x)≡−∂i(Aijxjq)+12Bij∂i∂jq (1)

and

 L†q(x)≡Aijxj∂iq+12Bij∂i∂jq (2)

where repeated upper and lower indices in lower case Latin are summed over. Thus is the component of the (linear) drift vector. We assume the diffusion matrix, whose elements are , is positive definite.

Throughout we assume that the left and right eigenvectors of form a complete set for , and that their corresponding eigenvalues (perhaps complex) have negative real part222Since is real, any complex eigenvectors occur as conjugate pairs.. These eigenvectors are normalized and bi-orthogonal

 wI⋅eJ=w∗I⋅e∗J=δIJ. (3)

### 2.1 Stationary States

The stationary state satisfies

 Lf0(x)=0.

and is given by

 f0(x)=1√(2π)NdetΣexp(−12xTΣ−1x). (4)

The covariance matrix is a solution to the Liapunov equation

 AΣ+ΣAT=−B. (5)

The corresponding stationary state of satisfies

 L†g0=0

and is given by

 g0(x)=1. (6)

As we will show, the eigenfunctions of and form a biorthogonal set, and for the stationary (ground) states we have

 ⟨g0,f0⟩≡∫g∗0(x)f0(x)dNx=1

as follows from the normalization of the multidimensional Gaussian .

### 2.2 Raising Operators

We generate the complete set of eigenfunctions of and by application of raising operators to and . The next two subsections address the two sets of eigenfunctions in turn.

#### 2.2.1 Forward Eigenfunctions

Let , be the right eigenvectors of with corresponding eigenvalues . (Recall we assume that all the eigenvalues of have negative real part, corresponding to an asympotitically-stable fixed point at .) Let denote the component of the eigenvector. The operators

 VI=−eiI∂i≡−eI⋅∇,I=1,2,…N (7)

satisfy the commutation relation

 [L,VI]=λIVI. (8)

(There is no sum on .) The commutator (8) establishes as a raising operator for eigenfunctions of : If is an eigenfunction of with eigenvalue , then is also an eigenfunction of with eigenvalue . The application of products of the to generate new eigenfunctions of which we denote by subscripts indicating the number of applications of each operator. For example

 fn1,n2,0,0…≡(V1)n1(V2)n2f0 with Lfn1,n2,0,0…=(n1λ1+n2λ2)fn1,n2,0,0…

Since , the subscript labels on uniquely determine the eigenfunctions. We build a general eigenfunction of as

 fn1,n2,…,nN≡(N∏I=1(VI)nI)f0 (9)

which satisfies

 (10)

Clearly from Eqn. (9)

 (VI)kfn1,…,nI,…,nN=fn1,…,(nI+k),…,nN. (11)

Let , be the left eigenvectors of with eigenvalues , and let denote its component. (The subscript indicating the component is appropriate since we regard as a co-vector.) Define the operators

 ¯VI≡2(w∗Iixi)+2[(A+λI)−1Bw∗I]i∂i=2[w∗Iixi−(Σw∗I)i∂i], (12)

where the second equality follows from the Liapunov equation (5). This operator satisfies the commutation relation

 [L†,¯VI]=λ∗I¯VI, (13)

which establishes it as a raising operator for eigenfunctions of : If is an eigenfunction of with eigenvalue , then is an eigenfunction with eigenvalue . Analogously to the forward eigenfunctions, acting with on generates a new eigenfunction which we denote by subscripts indicating the number of applications of each raising operator. This subscript notation is free of ambiguity since the raising operators corresponding to different eigenvectors of commute

 [¯VI,¯VJ]=−2(w∗J)TΣw∗I+2(w∗I)TΣw∗J=0.

We construct a general eigenfunction of by repeated application of raising operators. Thus

 gn1,n2,…,nN=N∏I=1(¯VI)nIg0, (14)

which satisfies

 L†gn1,n2,…,nN=(N∑I=1nIλ∗I)gn1,n2,…,nN. (15)

Clearly, from Eqn. (14)

 (¯VI)kgn1,…,nI,…,nn=gn1,…,(nI+k),…,nn. (16)

### 2.3 Lowering Operators

The two sets of raising operators are complemented by lowering operators. Their effect on the states, derived here, makes proving the bi-orthogonality property trivial.

Taking the adjoint of Eqn. (8) yields

 [L†,V†I]=−λ∗IV†I (17)

which establishes as a lowering operator for eigenfunctions of : If is an eigenfunction of with eigenvalue , then is an eigenfunction of with eigenvalue . In particular, kills the stationary state

 V†Ig0=e∗I⋅∇g0=0.

Using the commutators

 [V†I,¯VJ]=2e∗I⋅w∗J=2δIJ (18)

and , and the expression for the adjoint eigenfunctions Eqn. (14), the action of on eigenfunctions of follows as

 V†Jgn1,n2,…,nN ≡ V†J(N∏I=1(¯VI)nI)g0=⎛⎝∏I≠J(¯VI)nI⎞⎠V†J(¯VJ)nJg0 (19) = ⎛⎝∏I≠J(¯VI)nI⎞⎠[(¯VJ)nJV†J+(2nJ)(¯VJ)nJ−1]g0 = (2nJ)gn1,n2,…,nJ−1,…,nN

where the second line follows from the first by commuting past all factors of , and the third line follows since . Multiple applications yield

 (V†J)kgn1,n2,…,nN={2knJ!(nJ−k)!gn1,n2,…,nJ−k,…,nN,k≤nJ0,k>nJ. (20)

#### 2.3.2 Forward Eigenfunctions

Taking the adjoint of Eqn. (13) we recover

 [L,¯V†I]=−λI¯V†I (21)

which establishes as a lowering operator for eigenfunctions of : If is an eigenfunction of with eigenvalue , then is an eigenfunction with eigenvalue . In particular, it kills the stationary state

 ¯V†If0=2(wIixi+(ΣwI)i∂i)f0=0 (22)

having used the definition of in Eqn. (4). Similarly to the case for the adjoint lowering operators, it is straightforward to show that

 ¯V†Jfn1,n2,…,nN=(2nJ)fn1,…,nJ−1,…,nN. (23)

Multiple applications yield

 (¯V†J)kfn1,…,nN={cc2knJ!(nJ−k)!fn1,…,nJ−k,…,nN,k≤nJ0,k>nJ. (24)

## 3 Bi-Orthogonality and Normalization

The two sets of functions form a bi-orthogonal set. The usual result, that eigenfunctions corresponding to different eigenvalues are orthogonal can be strengthened. We will show that

 ⟨gm1,…mN,fn1,…,nN⟩ ≡ ∫g∗m1,…mN(x)fn1,…,nN(x)dNx (25) = (N∏I=1δmI,nI2nI(nI!)) = (N∏I=1δmI,nI2nI(nI!)).

We consider three cases:

##### Case I

Suppose for some . Write

 ⟨gm1,…,mJ,…,mN,fn1,…,nJ,…,nN⟩ = ⟨gm1,…,mJ,…,mN,VnJJfn1,…,nJ=0,…,nN⟩ (26) = ⟨(V†J)nJgm1,…,mJ,…,mN,fn1,…,nJ=0,…,nN⟩=0

having used Eqn. (20) to arrive at the last equality.

##### Case II

Now suppose instead that for some . Write

 ⟨gm1,…,mJ,…,mN,fn1,…,nJ,…,nN⟩ = ⟨¯VmJJgm1,…,mJ=0,…,mN,fn1,…,nJ,…,nN⟩ (27) = ⟨gm1,…,mJ=0,…,mN,(¯V†J)mJfn1,…,nJ,…,nN⟩=0

having used Eqn. (24) to arrive at the last equality.

##### Case III

Suppose for all . Write

 ⟨gn1,…,nN,fn1,…,nN⟩ = ⟨gn1,…,nN,VnJ1f0,n2,…,nN⟩ (28) = ⟨(V†1)n1gn1,n2,…,mN,f0,n2,…,nN⟩ = 2n1n1!⟨g0,n2,…,nN,f0,n2,…,nN⟩ = 2n1n1!⟨(V†2)n2g0,n2,n3,…,nN,f0,0,n3,…,nN⟩ = 2n12n2(n1!)(n2!)⟨g0,0,n3,…,nN,f0,0,n3,…,nN⟩ ⋮ = ⟨N∏J=12nJ(nJ!)⟩⟨g0,f0⟩ = ⟨N∏J=12nJ(nJ!)⟩

having used Eqn. (20) repeatedly. The three cases together prove the desired result Eqn. (25).

## 4 Analytic Form of the Eigenfunctions

In general coordinates, the eigenfunctions and do not assume a familiar analytic form. However they do if we transform to coordinates in which . (Since is a real, symmetric, positive-definite matrix, this is always possible.) In these coordinates,

 f0(x) = 1√πNe−xTx (29) g0(x) = 1. (30)

In these coordinates the forward eigenfunctions (9) are sums of products of Hermite polynomials times and the backward eigenfunctions (14) are sums of products of Hermite polynomials.

### 4.1 Forward Eigenfunctions

To start, note that the usual generating expression for the Hermite polynomials (Abramowitz and Stegun, 1972) can be rearranged to read

 Hn(xi)e−xTx=(−∂i)ne−xTx. (31)

Hence, application of to yields

 VIf0(x)=−eiI∂if0(x)=N∑i=1eiIH1(xi)f0(x)

a linear combination of Hermite functions in each of the variables . (When a coordinate index falls inside a function argument — as in the last expression — we will write the summation explicitly to avoid confusion.) Applying to , times results in

 VnIIf0=(−eiI∂i)nIf0

which can be evaluated using the multinomial theorem. Explicitly

 VnIIf0 = ∑(i1+…+iN=nI)(nI!i1!i2!⋯in!)∏1≤k≤N(−ekI∂k)ikf0 = ∑(i1+…+iN=iI)(nI!i1!i2!⋯in!)∏1≤k≤N(ekI)ikHik(xk)f0(x)

where the summation is over all values of indices satisfying the constraint .

The action of two distinct ladder operators multiple times is

 VnIIVnJJf0 = ∑\footnotesize(i1+…+iN=nI)(j1+…+jN=nJ)(nI!i1!i2!…iN!)(nJ!j1!j2!…jN!)∏1≤k≤N(ekI)ip(ekJ)jk(−∂k)ik+jkf0 = ∑\footnotesize(i1+…+iN=nI)(j1+…+jN=nJ)(nI!i1!i2!…iN!)(nJ!j1!j2!…jN!)∏1≤k≤N(ekI)ip(ekJ)jkHik+jk(xk)f0

This generalizes in the obvious way to products of the form

 fn1,n2,…,nN(x)=Vn11Vn22…VnNNf0(x). (33)

So in our special coordinates the eigenfunctions of are sums of products of Hermite polynomials times .

In our special coordinates so the backward raising operator (12) simplifies to

 ¯VI=2(w∗Iixi)−(w∗I)i∂i, (34)

where . From the standard recursion relations for the Hermite polynomials one has

 (2x−ddx)Hn(x)=Hn+1(x). (35)

Thus using the form of in our special coordinates (34), its action on a product of Hermite polynomials is

 ¯VIHi1(x1)Hi2(x2)⋯HiN(xN)=N∑k=1Hi1(x1)⋯w∗IkHik+1(xk)⋯HiN(xN).

It is convenient to define the operator which raises the order of by unity,

 rkHi1(x1)…Hik(xk)…HiN(xN)=Hi1(x1)…Hik+1(xk)…HiN(xN).

Then the action of on a product of Hermite polynomials is

 ¯VIHi1(x1)…HiN(xN) = N∑k=1w∗IkrkHi1(x1)…HiN(xN) (36) = N∑k=1Hi1(x1)…w∗IkHik+1(xk)…HiN(xN).

Writing

 g0=1=H0(x1)H0(x2)⋯H0(xN)

and using Eqn. (36) gives

 ¯VIg0=N∑k=1w∗IkH1(xk) (37)

The action of on is evaluated using the multinomial theorem

 ¯VnIIg0 = ∑\footnotesize(i1+⋯+in=nI)(nI!i1!…iN!)∏1≤k≤N(w∗Ikrk)ikg0 = ∑\footnotesize(i1+⋯+in=nI)(nI!i1!…iN!)∏1≤k≤N(w∗Ik)ikHik(xk)

The action of two such operators on is

 ¯VnJJ¯VnIIg0 = ∑\footnotesize(i1+⋯+iN=nI)(j1+⋯+jN=nJ)(nJ!j1!…jN!)(nI!i1!…iN!)∏1≤k≤N(w∗Ik)ik(w∗Jk)jkHik+jk(xk).

This generalizes in the obvious way to evaluate the general backward eigenfunction in Eqn. (14).

## 5 Examples

##### Example 1 — Derivative and Multiplication Operators

Recall that the left () and right () eigenvectors of form a complete, bi-orthogonal set for (see Eqn. (3)). Then, from the definition of the raising operator (7) for eigenfunctions of , we have

 ∇=N∑I=1w∗IV†I. (39)

Similarly, from the definition of the raising operator (12) for eigenfunctions of and this last result (39), we have

 x=12N∑I=1e∗I[¯VI+2N∑J=1(w∗I⋅Σw∗J)V†J]. (40)
##### Example 2 — Forward and Adjoint FPE Operators

Using the results in Example 1, we can rewrite as

 L=12∑IλIVI¯V†I (41)

which can be verified by its action on the eigenfunctions

 Lfn1,…,nN = 12∑IλIVI¯V†Ifn1,…,nN=∑InIλIVIfn1,…,nI−1,…,nN (42) = ∑InIλIfn1,…,nN.

Similarly, we recover as

 L†=12∑Iλ∗I¯VIV†I (43)

which can be verified by its action on the eigenfunctions . (Both Eqns. (41) and (43) have the flavor of quantum oscillator Hamiltonians involving products of lowering and raising operators333Hence, and are occupation number operators for the eigenstate in and respectively. The difference here is that and are not Hermitian. The correspondence with quantum mechanics (in the 1-D case) is discussed by Gardiner (2009).)

##### Example 3 — Fourier Series

Time-dependent solutions can be expanded in terms of the forward eigenfunctions exactly as in the one-dimensional case. Let denote an index set for the eigenfunctions and eigenvalues. Then

 F(x,t)=∑KαKfK(x)exp(λKt) (44)

clearly satisfies

 ∂tF(x,t)=LF(x,t)

where the coefficients are given by . (Since the eigenvalues of have negative real part, all the exponentials in (44) are decaying.) In multiple dimensions, the drift Jacobian can have one or more pairs of complex-conjugate eigenvalues and Eqn. (44) can represent damped oscillating solutions.

##### Example 4 — Solutions to Inhomogeneous Equations

In our perturbation solution for densities satisfying a nonlinear master equation (Leen and Friel, 2011, 2012; Leen et al., 2012), the order corrections to the equilibrium density is denoted and is given by inhomogeneous equations of the form

 LP(i)(x)=q(i)(x), (45)

where the are known. Let denote an index set for the eigenfunctions and eigenvalues. Next, expand in a linear combination of the eigenfunctions , substitute that into (45), take the inner product with , and use bi-orthogonality relations (25) to obtain

 P(i)(x)=∑K≠0⟨gK,q(i)⟩λK⟨gK,fK⟩fK(x) (46)

where the summation excludes .

## 6 Conclusion

We have provided raising and lowering operators to develop the eigenfunctions (and their corresponding eigenvalues) of the forward and adjoint multidimensional Fokker-Planck operators for the Ornstein-Uhlenbeck process. The eigenfunctions form a basis for expanding solutions to the time-dependent Fokker-Planck equation, and for a perturbation expansion of the densities arising from a nonlinear master equation (Leen and Friel, 2011, 2012; Leen et al., 2012; Thomas and Grima, 2015).

We gave bi-orthogonality and normalization results. We showed that in coordinates for which the covariance of the stationary state is spherically symmetric with variance one-half, the eigenfunctions reduce to sums of products of Hermite polynomials times . In applications to time-dependent solutions of the Fokker-Planck equation and to inhomogeneous equations (see Example 3 and Example 4 in Section 5) one assumes the eigenfunctions form a complete set on . The proof of completeness is similar to that used to show completeness of the one-dimensional Hermite functions.

##### Acknowledgments

This work was supported by NSF under grant IIS-0812687. The authors thank Gerardo Lafferriere and Crispin Gardiner for their comments.

## References

• Abramowitz and Stegun (1972) Abramowitz, M., Stegun, L., 1972. Handbook of Mathematical Functions. U.S. Department of Commerce, National Bureau of Standards.
• Gardiner (2009) Gardiner, C., 2009. Stochastic Methods, A Handbook for the Natural and Social Sciences, Fourth Edition. Springer-Verlag, Berlin.
• Leen and Friel (2011) Leen, T.K., Friel, R., 2011. Perturbation theory for stochastic learning dynamics, in: Proceedings of the IJCNN 2011, IEEE Press, San Jose, CA.
• Leen and Friel (2012) Leen, T.K., Friel, R., 2012. Stochastic perturbation methods for spike-timing-dependent plasticity. Neural Computation. 24, 1109–1146.
• Leen et al. (2012) Leen, T.K., Friel, R., Nielsen, D., 2012. Approximating distributions in stochastic learning. Neural Networks 32, 219–228.
• Liberzon and Brockett (2000) Liberzon, D., Brockett, R.W., 2000. Spectral analysis of Fokker-Planck and related operators arising from linear stochastic differential equation. SIAM J. Control Optim. 38, 1453–1467.
• Risken (1989) Risken, H., 1989. The Fokker-Planck Equation. Springer-Verlag, Berlin.
• Thomas and Grima (2015) Thomas, P., Grima, R., 2015. Approximate probability distributions of the master equation. Physical Review E92. DOI: 10.1103/PhysRevE.92.012120.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters