Processing math: 100%
Press "Enter" to skip to content

Random walk, Dirichlet problem, and Gaussian free field

Solution of a Dirichlet problem obtained with a Monte Carlo simulation, using the Julia programming language
Solution of a Dirichlet problem by Monte Carlo simulation of probabilistic formulation, using Julia.

This post is about the discrete Dirichlet problem and Gaussian free field, linked with the random walk on Zd. Discreteness allows to go to the concepts with minimal abstraction. The Gaussian free field (GFF) is an important Gaussian object which appears, like most Gaussian objects, as a limiting object in many models of Statistical Physics.

Symmetric Nearest Neighbors random walk. The symmetric nearest neighbors random walk on Zd is the sequence of random variables X=(Xn)n0 defined by the linear recursion

Xn+1=Xn+ξn+1=X0+ξ1++ξn+1

where (ξ)n1 is a sequence of independent and identically distributed random variables on Zd, independent of X0, and uniformly distributed on the discrete 1 sphere: {±e1,,±ed} where e1,,ed is the canonical basis of Rd. The term symmetric comes from the fact that in each of the d dimensions, both directions are equally probable in the law of the increments.

The sequence X is a homogeneous Markov chain with state space Zd since

P(Xn+1=xn+1|Xn=xn,,X0=x0)=P(Xn+1=xn+1|Xn=xn)=P(X1=xn+1|X0=xn)=P(ξ1=xn+1xn).

Its transition kernel P:Zd×Zd[0,1] is given for every x,yZd by

P(x,y)=P(Xn+1=y|Xn=x)=1{|xy|1=1}2dwhere|xy|1=dk=1|xkyk|.

It is an infinite Markov transition matrix. It acts on a bounded function f:ZdR as

(Pf)(x)=yZdP(x,y)f(y)=E(f(X1)|X0=x).

The sequence X is also a martingale for the filtration (Fn)n0 defined by F0=σ(X0) and Fn+1=σ(X0,ξ1,,ξn+1). Indeed, by measurability and independence,

E(Xn+1|Fn)=Xn+E(ξn+1|Fn)=Xn+E(ξn+1)=Xn.

Dirichlet problem. For any bounded function f:ZdR and any xZd,

E(f(Xn+1)|Xn=x)f(x)=(Pf)(x)f(x)=(Δf)(x)whereΔ=PI.

Here I is the identity operator I(x,y)=1{x=y}. The operator

Δ=PI

is the generator of the symmetric nearest neighbors random walk. It is a discrete Laplace operator, which computes the difference between the value at a point and the mean on its neighbors for the 1 distance:

(Δf)(x)={12dy:|yx|1=1f(y)}f(x)=12dy:|yx|1=1(f(y)f(x)).

We say that f is harmonic on AZd when Δf=0 on A, which means that on every point of A, the value of f is equal to the mean of the values of f on the 2d neighbors of x. The operator Δ is local in the sense that (Δf)(x) depends only on the value of f at x and its nearest neighbors. Consequently, the value of Δf on A depends only on the values of f on

ˉA:=AA

where

A:={yA:xA,|xy|1=1}

is the exterior boundary of A. In this context, the Dirichlet problem consists in finding a function which is harmonic on A and for which the value on A is prescribed. This problem is actually a linear algebra problem for which the result below provides a probabilistic expression of the solution, based on the stopping time

τA:=inf{n0:XnA}.

In Physics, the quantity f(x) models typically the temperature at location x, while the harmonicity of f on A and the boundary condition on A express the thermal equilibrium and a thermostatic boundary, respectively.

Theorem 1 (Dirichlet problem) Let AZd be finite. Then for any xA,

Px(τA<)=1.

Moreover, for any g:AR, the function f:ˉAR defined for any xˉA by

f(x)=Ex(g(XτA))

is the unique solution of the system

{f=gon A,Δf=0on A.

When d=1 we recover the function which appears in the gambler ruin problem. Recall that the image of a Markov chain by a function which is harmonic for the generator of the chain is a martingale (discrete Itô formula !), and this explains a posteriori the probabilistic formula for the solution of the Dirichlet problem, thanks to the Doob optional stopping theorem.

The quantity f(x)=Ex(g(τA))=yAg(y)Px(τA=y) is the mean of g for the law μx on A called the harmonic measure defined for any yA by

μx(y):=Px(XτA=y)

Proof: The property Px(τA<)=1 follows from the Central Limit Theorem or by using conditioning which gives a geometric upper bound for the tail of τA.

Let us check that the proposed function f is a solution. For any xA, we have τA=0 on {X0=x} and thus f=g on A. Let us show now that Δf=0 on A. We first reduce the problem by linearity to the case g=1{z} with zA. Next, we write, for any yˉA,

f(y)=Py(XτA=z)=n=0Py(Xn=z,τA=n)=1y=z+1yAn=1x1,,xn1AP(y,x1)P(x1,x2)P(xn1,z).

On the other hand, since f=0 on A{z} and since Δ is local, we have, for any xA,

(Pf)(x)=yZdP(x,y)f(y)=P(x,z)f(z)+yAP(x,y)f(y).

As a consequence, for any xA,

(Pf)(x)=P(x,z)f(z)+n=1y,x1,,xn1AP(x,y)P(y,x1)P(xn1,z)=P(x,z)f(z)+(f(x)(1x=z+P(x,z))),

which is equal to f(x) since 1x=z=0 and f(z)=1. Hence Pf=f on A, in other words Δf=0 on A.

To establish the uniqueness of the solution, we first reduce the problem by linearity to show that f=0 is the unique solution when g=0. Next, if f:ˉAR is harmonic on A, the interpretation of Δf(x) as a difference with the mean on the nearest neighbors allows to show that both the minimum and the maximum of f on ˉA are (at least) necessarily achieved on the boundary A. But since f is vanishes on A, it follows that f vanishes on ˉA. ☐

Dirichlet problem and Green function. Here is a generalization of Theorem 1. When h=0 we recover Theorem 1.

Theorem 2 (Dirichlet problem and Green function) If AZd is finite then for any g:AR and h:AR, the function f:ˉAR defined for any xˉA by

f(x)=Ex(g(XτA)+τA1n=0h(Xn))

is the unique solution of

{f=gon A,Δf=hon A.

When d=1 and h=1 then we recover a function which appears in the analysis of the gambler ruin problem.

For any xA we have

f(x)=yAg(y)Px(τA=y)+yAh(y)GA(x,y)

where GA(x,y) is the average number of passages at y starting from x and before escaping from A:

GA(x,y):=Ex{τA1n=01Xn=y}=n=0Px(Xn=y,n<τA).

We say that GA is the Green function of the symmetric nearest neighbors random walk on A killed at the boundary A. It is the inverse of the restriction ΔA of Δ to functions on ˉA vanishing on A:

GA=Δ1A.

Indeed, when g=0 and h=1{y} we get f(x)=GA(x,y) and thus ΔAGA=IA.

Proof: Thanks to Theorem 1 it suffices by linearity to check that f(x)=1xAGA(x,z) is a solution when g=0 and h=1{z} with zA. Now for any xA,

f(x)=1{x=z}+n=1Px(Xn=z,n<τA)=1{x=z}+y:|xy|1=1n=1P(Xn=z,n<τA|X1=y)P(x,y)=1{x=z}+u:|xy|1=1f(y)P(x,y)

thanks to the Markov property. We have f=h+Pf on A, in other words Δf=h on A. ☐

Beyond the symmetric nearest neighbors random walk. The proofs of Theorem 1 and Theorem 2 remain valid for asymmetric nearest neighbors random walks on Zd, provided that we replace the generator Δ of the symmetric nearest neighbors random walk by the generator L:=PI, which is still a local operator: |xy|>1L(x,y)=P(x,y)=0. One may even go beyond this framework by adapting the notion of boundary: {yA:xA,L(x,y)>0}.

Gaussian free field. The Gaussian free field is model of Gaussian random interface for which the covariance matrix is the Green function of the discrete Laplace operator. More precisely, let AZd be a finite set. An interface is a height function f:ˉAR which associates to each site xˉA a height f(x), also called spin in the context of Statistical Physics. For simplicity, we impose a zero boundary condition: f=0 on the boundary A of A.

Let FA be the set of interfaces f on ˉA vanishing at the boundary A, which can be identified with RA. The energy HA(f) of the interface fFA is defined by

HA(f)=14d{x,y}ˉA|xy|1=1(f(x)f(y))2,

where f=0 on A. The flatter is the interface, the smaller is the energy HA(f). Denoting u,vA:=xAu(x)v(x), we get, for any fFA,

HA(f)=14dxAyˉA|xy|1=1(f(x)f(y))f(x)=12Δf,fA.

Since HA(0)=0 and HA(f)=0 give f=0, the quadratic form HA is not singular and we can define the Gaussian law QA on FA which favors low energies:

QA(df)=1ZAeHA(f)dfwhereZA:=FAeHA(f)df.

This Gaussian law, called the Gaussian free field, is characterized by its mean mA:AR and its covariance matrix CA:A×AR, given for any x,yA by

mA(x):=fxQA(df)=0

and

CA(x,y):=fxfyQA(df)mA(x)mA(y)=Δ1A(x,y)=GA(x,y),

where fx denotes the coordinate map fx:fFAf(x)R.

Simulation of the Gaussian free field. The simulation of the Gaussian free field can be done provided that we know how to compute a square root of GA in the sense of quadratic forms, using for instance a Cholesky factorization or diagonalization. Let us consider the square case:

A=L(0,1)dZd={1,,L1}d.

This gives ˉA={0,1,,L}d. In this case, one can compute the eigenvectors and eigenvalues of ΔA and deduce the ones of GA=Δ1A. In fact, the continuous Laplace operator Δ on [0,1]dRn, with Dirichlet boundary conditions, when defined on the Sobolev space H20([0,1]d), is a symmetric operator with eigenvectors {en:n{1,2,}d} and eigenvalues {λn:n{1,2,}d} given by

en(t):=2d/2di=1sin(πniti),andλn:=π2|n|22=π2(n21++n2d).

Similarly, the discrete Laplace operator ΔA with Dirichlet boundary conditions on this A has eigenvectors {eLn:nA} and eigenvalues {λLn:nA} given for any kA by

eLn(k):=en(kL)andλLn:=2di=1sin2(πni2L).

The vectors eLn vanish on the boundary A. It follows that for any x,yA,

ΔA(x,y)=nAλLneLn,fxeLn,fy=nAλLneLn(x)eLn(y)

and

GA(x,y)=nA(λLn)1eLn(x)eLn(y).

A possible matrix square root GA of GA is given for any x,yA by

GA(x,y)=nA(λLn)1/2eLn(x)eLn(y).

Now if Z=(Zy)yA are independent and identically distributed Gaussian random variables with zero mean and unit variance, seen as a random vector of RA, then the random vector

GAZ=(yAGA(x,y)Zy)xA=(yAnA(λLn)1/2eLn(x)eLn(y)Zy)xA

is distributed according to the Gaussian free field QA.

Continuous analogue. The scaling limit εZdRd leads to natural continuous objects (in time and space): the simple random walk becomes the standard Brownian Motion via the Central Limit Theorem, the discrete Laplace operator becomes the Laplace second order differential operator via a Taylor formula, while the discrete Gaussian free field becomes the continuous Gaussian free field (its covariance is the Green function of the Laplace operator). The continuous object are less elementary since their require much more analytic and probabilistic machinery, but in the mean time, they provide differential calculus, a tool which is problematic in discrete settings due in particular to the lack of chain rule.

Further reading.

Note. The probabilistic approach to the Dirichlet problem goes back at least to Shizuo Kakutani (1911 – 2004), who studied around 1944 the continuous version with Brownian Motion. This post is the English translation of an excerpt from a forthcoming book on stochastic models, written, in French, together with my old friend Florent Malrieu.

1D GFF simulated with Julia
1D GFF simulated with Julia
2D GFF simulated with Julia
2D GFF simulated with Julia
Johann Peter Gustav Lejeune Dirichlet (1805 - 1859) who probably never thought about random walks
Johann Peter Gustav Lejeune Dirichlet (1805 – 1859) who probably never thought about random walks and Gaussian fields
#
# load within julia using include("dirichlet.jl")
#
nx, ny = 25, 25
f = zeros(nx,ny)
f[1,:], f[nx,:], f[:,1], f[:,ny] = 2, 2, -2, -2
C = zeros(nx,ny)
N = 5000
#
for x = 1:nx
 for y = 1:ny
  for k = 1:N
   xx = x
   yy = y
   while ((xx > 1) && (xx < nx) && (yy > 1) && (yy < ny))
    b = rand(2)
    if (b[1] < .5)
     if (b[2] < .5)
      xx -= 1
     else
      xx += 1
     end
    elseif (b[2] < .5)
     yy -= 1
    else
     yy += 1
    end
   end # while
   C[x,y] += f[xx,yy]
  end # for k
 end # for y
end # for x
using Compose, Gadfly
draw(PNG("dirichlet.png", 10cm, 10cm), spy(C/N))
#
# load within julia using include("gff.jl")
#
Pkg.add("Gadfly") # produces "INFO: Nothing to be done." if already installed.
Pkg.add("Cairo") # idem.
#Pkg.add("PyPlot") # idem.
using Compose, Gadfly, Cairo # for 2D graphics
#using PyPlot # for 3D graphics (Gadfly does not provide 3D graphics)
#
# dimension 1
#
L = 500
Z = randn(L-1)
lambda = 2*sin(pi*[1:L-1]/(2*L)).^2
evec = sqrt(2)*sin(pi*kron([1:L-1],[1:L-1]')/L)
GFF1 = zeros(L-1)
for x = 1:L-1
    for y = 1:L-1
        for n = 1:L-1
            GFF1[x] += (lambda[n])^(-1/2)*evec[n,x]*evec[n,y]*Z[y]
        end
    end
end
# graphics ith Gadfly
p = plot(x = [1:L-1], y = GFF1,
         Geom.line,
         Guide.xlabel(""),
         Guide.ylabel(""),
         #Guide.title("GFF"),
         Theme(default_color=color("black")))
draw(PNG("gff1.png", 10cm, 10cm), p)
#
# dimension 2
#
L = 20
Z = randn(L-1,L-1)
AUX = kron(ones(L-1)',[1:L-1])
lambda = 2*sin(pi*AUX/(2*L)).^2 + 2*sin(pi*AUX'/(2*L)).^2
evec = zeros(L-1,L-1,L-1,L-1) # n1,n2,k1,k2 : e_n^L(k)
for n1 = 1:L-1
 for n2 = 1:L-1
  AUX = kron(ones(L-1)',[1:L-1])
  evec[n1,n2,:,:] = reshape(2*sin(pi*n1*AUX/L).*sin(pi*n2*AUX'/L),(1,1,L-1,L-1))
 end
end
GFF2 = zeros(L-1,L-1)
for x1 = 1:L-1
 for x2 = 1:L-1
  for y1 = 1:L-1
   for y2 = 1:L-1
    for n1 = 1:L-1
     for n2 = 1:L-1
      aux = (lambda[n1,n2])^(-1/2)
      aux *= evec[n1,n2,x1,x2]
      aux *= evec[n1,n2,y1,y2]
      aux *= Z[y1,y2]
      GFF2[x1,x2] += aux
     end
    end
   end
  end
 end
end
# graphics with PyPlot
#surf(GFF2)
#savefig("gff2.svg")
# graphics with Gadfly
draw(PNG("gff2.png", 10cm, 10cm), spy(GFF2))
    Leave a Reply

    Your email address will not be published.

    This site uses Akismet to reduce spam. Learn how your comment data is processed.

    Syntax · Style · .