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)n≥0 defined by the linear recursion
Xn+1=Xn+ξn+1=X0+ξ1+⋯+ξn+1
where (ξ)n≥1 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+1−xn).
Its transition kernel P:Zd×Zd→[0,1] is given for every x,y∈Zd by
P(x,y)=P(Xn+1=y|Xn=x)=1{|x−y|1=1}2dwhere|x−y|1=d∑k=1|xk−yk|.
It is an infinite Markov transition matrix. It acts on a bounded function f:Zd→R as
(Pf)(x)=∑y∈ZdP(x,y)f(y)=E(f(X1)|X0=x).
The sequence X is also a martingale for the filtration (Fn)n≥0 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:Zd→R and any x∈Zd,
E(f(Xn+1)|Xn=x)−f(x)=(Pf)(x)−f(x)=(Δf)(x)whereΔ=P−I.
Here I is the identity operator I(x,y)=1{x=y}. The operator
Δ=P−I
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)={12d∑y:|y−x|1=1f(y)}−f(x)=12d∑y:|y−x|1=1(f(y)−f(x)).
We say that f is harmonic on A⊂Zd 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:=A∪∂A
where
∂A:={y∉A:∃x∈A,|x−y|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{n≥0:Xn∈∂A}.
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 ∅≠A⊂Zd be finite. Then for any x∈A,
Px(τ∂A<∞)=1.
Moreover, for any g:∂A→R, the function f:ˉA→R 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))=∑y∈∂Ag(y)Px(τ∂A=y) is the mean of g for the law μx on ∂A called the harmonic measure defined for any y∈∂A 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 x∈∂A, 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 z∈∂A. Next, we write, for any y∈ˉA,
f(y)=Py(Xτ∂A=z)=∑∞n=0Py(Xn=z,τ∂A=n)=1y=z+1y∈A∑∞n=1∑x1,…,xn−1∈AP(y,x1)P(x1,x2)⋯P(xn−1,z).
On the other hand, since f=0 on ∂A∖{z} and since Δ is local, we have, for any x∈A,
(Pf)(x)=∑y∈ZdP(x,y)f(y)=P(x,z)f(z)+∑y∈AP(x,y)f(y).
As a consequence, for any x∈A,
(Pf)(x)=P(x,z)f(z)+∑∞n=1∑y,x1,…,xn−1∈AP(x,y)P(y,x1)⋯P(xn−1,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:ˉA→R 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 ∅≠A⊂Zd is finite then for any g:∂A→R and h:A→R, the function f:ˉA→R defined for any x∈ˉA by
f(x)=Ex(g(Xτ∂A)+τ∂A−1∑n=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 x∈A we have
f(x)=∑y∈∂Ag(y)Px(τ∂A=y)+∑y∈Ah(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{τ∂A−1∑n=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)=1x∈AGA(x,z) is a solution when g=0 and h=1{z} with z∈A. Now for any x∈A,
f(x)=1{x=z}+∑∞n=1Px(Xn=z,n<τ∂A)=1{x=z}+∑y:|x−y|1=1∑∞n=1P(Xn=z,n<τ∂A|X1=y)P(x,y)=1{x=z}+∑u:|x−y|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:=P−I, which is still a local operator: |x−y|>1⇒L(x,y)=P(x,y)=0. One may even go beyond this framework by adapting the notion of boundary: {y∉A:∃x∈A,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 ∅≠A⊂Zd be a finite set. An interface is a height function f:ˉA→R 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 f∈FA is defined by
HA(f)=14d∑{x,y}⊂ˉA|x−y|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,v⟩A:=∑x∈Au(x)v(x), we get, for any f∈FA,
HA(f)=14d∑x∈A∑y∈ˉA|x−y|1=1(f(x)−f(y))f(x)=12⟨−Δf,f⟩A.
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)=1ZAe−HA(f)dfwhereZA:=∫FAe−HA(f)df.
This Gaussian law, called the Gaussian free field, is characterized by its mean mA:A→R and its covariance matrix CA:A×A→R, given for any x,y∈A 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:f∈FA↦f(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)d∩Zd={1,…,L−1}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]d∈Rn, 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/2d∏i=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:n∈A} and eigenvalues {λLn:n∈A} given for any k∈A by
eLn(k):=en(kL)andλLn:=−2d∑i=1sin2(πni2L).
The vectors eLn vanish on the boundary ∂A. It follows that for any x,y∈A,
ΔA(x,y)=∑n∈AλLn⟨eLn,fx⟩⟨eLn,fy⟩=∑n∈AλLneLn(x)eLn(y)
and
GA(x,y)=−∑n∈A(λLn)−1eLn(x)eLn(y).
A possible matrix square root √GA of GA is given for any x,y∈A by
√GA(x,y)=−∑n∈A(λLn)−1/2eLn(x)eLn(y).
Now if Z=(Zy)y∈A 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=(∑y∈A√GA(x,y)Zy)x∈A=(∑y∈A∑n∈A(−λLn)−1/2eLn(x)eLn(y)Zy)x∈A
is distributed according to the Gaussian free field QA.
Continuous analogue. The scaling limit εZd→Rd 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.
- James Glimm and Arthur Jaffe. Quantum Physics: A Functional Integral Point of View, mostly chapter 7 “Covariance Operator = Green’s Function = Resolvent Kernel = Euclidean Propagator = Fundamental Solution” (great source of inspiration);
- Richard F. Bass. Probabilistic techniques in analysis (very pleasant)
- Scott Scheffield. Gaussian free fields for mathematicians (not so bad, after all)
- Gregory Lawler. Intersections of random walks mostly chapter 1 (old memories from my Master degree memoir!).
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.
# # 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))