Nonlinear Dimensionality Reduction II: Diffusion Maps

1 .Nonlinear Dimensionality Reduction II: Diffusion Maps 36-350, Data Mining 9 October 2009 Contents 1 Introduction 1 2 Diffusion-Map Coordinates 3 2.1 Fun with Transition Matrices . . . . . . . . . . . . . . . . . . . . 5 2.2 Multiple Scales . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Choosing q . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3 What to Do with the Diffusion Map Once You Have It 9 3.1 Spectral Clustering . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.2 Asymmetry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4 The Kernel Trick 11 A From Random Walks to the Diffusion Equation 14 1 Introduction Let’s re-cap what we did with locally linear embedding. We start p-dimensional data (x1 , . . . xn in vector form, stacked into a data-frame X), which we suspect are or are near a q-dimensional manifold embedded in the feature space. As an exercise in least squares, we got weights wij saying how much of data-point xj we needed to use to reconstruct xi . Since j wij = 1, we can think of wij as saying how similar i is to j. Next, we asked for the new vectors yi which were best reconstructed by those weights, i.e., which minimized 2 n n yi − wij yj (1) i=1 j=1 under the constraint that n−1 YT Y = I. We turned this into an eigenvalue problem, that of finding the eigenvectors of M = (I − w)T (I − w) (2) 1 

2 .It was nice to do this, because finding eigenvectors is something our computers are very good at. These eigenvectors were then the new coordinates of our data points on the low-dimensional manifold. To introduce diffusion maps, we need to change tack a little from the ge- ometric/optimization approach we’ve been taking, though we’ll come back to that. Notice that when we do LLE, we identify a neighborhood of k points for each point in our data set. We can imagine drawing a graph or network which connects a point to its neighbors. We can also attach weights to the edges in this graph, since we get them from our w matrix. Now, we’ve seen networks with weighted edges before in this class, when we looked at page-rank: there the edge-weights reflected how many links there were between pages, and they were normalized to sum to one for each page. It is not at all a coincidence that each row in w also sums to one. If the weights are all non-negative, we can use them to define a random walk on the data, just as in page-rank we had a random walk on the web.1 Let’s be a little more formal, and more general. Suppose that we have a function K which takes two data-points x1 , x2 are inputs, and gives us a non- negative real number as an output: K(x1 , x2 ) ≥ 0 (3) Further suppose that it’s symmetric K(x2 , x1 ) = K(x1 , x2 ) (4) We’ll form a matrix K for our data, where Kij = K(xi , xj ) (5) and as a last requirement, insist that the matrix be non-negative, that for any vector v v T Kv ≥ 0 (6) Tradition says that we should call this the kernel matrix or Gram matrix corresponding to the kernel2 K. The point of this matrix is, once again, to give us a weighted graph, which tells us which points are tied to which other points, and how similar those neighbors are — it’s going to play much the same role in the generalized case that w did for LLE.3 The advantage of abstracting away from the least-squares procedure of LLE to “some kernel function or other” is that it lets us deal with cases where least-squares does not apply — because it’s not the right metric, because some of the data are categorical, etc.; I will return to the “kernel trick” below, and we’ll come back to it again later in the course. 1 We didn’t require the weights to be non-negative, but they often will be spontaneously, and it’s not, generally, a hard requirement to impose. See Saul and Roweis (2003), if you want the details. 2 A name whose origins are lost in the mists of 19th century algebra. 3 If you are getting suspicious because w is not symmetric, even if we force is entries to be non-negative, good; but notice that what ultimately shows up in finding the coordinates is wT + w + wT w, which is symmetric. 2 

3 . Given this K, we need one more step to make a random walk on the data graph, which is to turn it into a stochastic transition matrix, where each row adds up to one: D = diag(rowSums(K)) (7) that is, D is the diagonal matrix whose ith entry is the sum of the ith row of K.4 Now D−1 K (8) defines the same graph as K, but it’s properly row-normalized. Let’s call this matrix A. It’s the transition matrix for a Markov chain — specifically or a random walk on the graph of the data. The random walk is biased towards making transitions that take it between similar points. If we look at the behavior of the random walk, then, it should tell us about the geometry of the data (how can we move around the data without ever making a big, discontinuous step?), and about groups of similar points (i.e., what are the clusters?). We’ll take these two in order. 2 Diffusion-Map Coordinates Let’s say that we want to assign a single coordinate to every point on the data, which we’ll write as f , with f (i) being the coordinate for i. Because there are n data points, instead of treating f as a function, we can also treat it as an n × 1 matrix, which for simplicity I’ll also write f . (Here, as with PCA and with LLE, starting with a single coordinate, the q = 1 case, is just to simplify the initial algebra.) Let’s say that we want to minimize Φ(f ) ≡ Aij (fi − fj )2 (9) i,j This should force the coordinates of points which are very similar (Aij ≈ 1) to be close to each other, but let the coordinates of deeply dissimilar points (Aij ≈ 0) vary freely.5 Now that we know what to do with quadratic forms, let’s try and turn this into one. Aij (fi − fj )2 = Aij (fi2 − 2fi fj + fj2 ) (10) i,j ij = fi2 Aij − 2fi Aij fj + Aij fj2 (11) i j j j = fi2 − 2 fi Aij fj + fj2 (12) i i j j = 2f T f − 2f T Af (13) 4 If you don’t like my mix of R and matrix notation, you try writing it out in proper matrix form. 5 This is not quite the same as the minimization problem for LLE. 3 

4 .using the fact that Aij = Aji . So minimizing Φ is the same as minimizing f T (I − A)f ≡ f T Lf (14) where of course L ≡ I − A. The matrix L is called the Laplacian of the graph. Now we impose the constraint that f T f = 1, to avoid the uninteresting minimum at f = 0, and as before we get an eigenvalue problem: Lf = λf (15) is the solution. Substituting back in to the expression for Φ, we see that the value of Φ is 2λ, so we want the bottom eigenvectors, as we did with LLE. For exactly that same reasons that w1 = 1 for LLE, we know that A1 = 1 — i.e., the vector of all 1s is an eigenvector of A with eigenvalue 1. It turns out (see below) that this is the largest eigenvalue. Now, every eigenvector of A is also an eigenvector of L, but the eigenvalues change: Lf = λf (16) (I − A)f = λf (17) f − Af = λf (18) Af = (1 − λ)f = µf (19) where of course µ = 1 − λ. Thus, 1 is an eigenvector of the Laplacian with eigenvalue 0. Again as with LLE, we discard this solution as useless. The next-to-bottom eigenvector of L has to be orthogonal to 1, because all the eigenvectors are orthogonal. This means it must have both negative and positive entries (otherwise 1T f = i fi > 0). We will see later how to use the signs of the entries in this eigenvector. What matters for right now is that it gives us a non-trivial coordinate, in which similar data points are close to each other. If we want q coordinates, then (yet again as with LLE) we just take the bottom q +1 eigenvectors6 f (n) , f (n−1) , . . . f (n−q) of L, and their eigenvalues 0 = λn , λn−1 , . . . λn−q . We discard the bottom one as uninteresting. The diffusion map Ψ uses the remaining eigenvectors and eigenvalues to find a point in Rq corresponding to each of the original data points. Specifically, the image of the point v is Ψ(v) = µ1 fv(n−1) , µ2 fv(n−2) , . . . µq f (n−q) (20) where I’ve abbreviated 1 − λn−i as µi . The coordinates we get from the diffusion map are related to, but not iden- tical with, the coordinates we would get from LLE.7 They former are, however, in a reasonable sense the optimal coordinates to represent the original matrix of similarities K (Lee and Wasserman, 2008). By now you should be asking what any of this has to do with “diffusion”. 6 Eigenvectors are conventionally numbered starting from 1 at the top. 7 In fact, in some cases, it can be shown (Belkin and Niyogi, 2003, §5) that the matrix in the LLE minimization problem is related to the Laplacian, because (I − w)T (I − w) ≈ 21 L2 . Since the powers of L have the same eigenvectors as L, when this holds the coordinates we get from the diffusion map are approximately the same as the LLE coordinates. 4 

5 .2.1 Fun with Transition Matrices The matrix A is a stochastic transition matrix, so the entries in each row are non-negative and add up to one. Suppose we have a function f on the graph; since there are only n data points, we can represent it as an n × 1 matrix, which for simplicity I’ll also write as f . What’s Af ? Well, it’s another n × 1 matrix: n (Af )i = Aij fj (21) j=1 In words, the ith entry of Af is a weighted average of the values of f at i’s neighbors; the weight of j is the likelihood of going from i to j. Suppose I write Zt for the node occupied by the random walk at time t. Then E [f (Zt+1 )|Zt = i] = (Af )i (22) If f encodes a function, then Af encodes the expected value of that function after one step of the random walk, A2 f is the expected function after two steps, and so forth. So, acting to the right, A forecasts what the value of a fixed function is going to be later. This is inside-out from (or, a the mathematicians say, “adjoint to”) asking where the random walk will be later. As you know, if ρ is a distribution over the states of a Markov chain, written as a 1 × n matrix, then ρA (23) is the distribution of the Markov chain at the next time step. So from the left A updates distributions, and from the right A forecasts functions. Without getting in to the more profound aspects of this duality8 , what we care about here is how this can help us understand the data. We know (if only because I asserted it during the lecture on page-rank) that A has a left eigenvector with eigenvalue 1: ρ(1) A = ρ(1) (24) Moreover, every entry of ρ(1) is > 0, and it is the top (or “dominant”) eigen- vector, the one with the largest eigenvalue. This means that all the other eigenvectors must have both positive and negative entries. What happens if we let the walk evolve for many steps undisturbed? Er- godicity happens, that’s what. Let’s say that the let eigenvectors of A are ρ(1) , ρ(2) , . . . ρ(n) , with corresponding eigenvalues 1 = µ1 > µ2 > . . . µn > 0. Since th eigenvectors are orthogonal, and there are n of them, they span the space — any arbitrary vector can be written as a linear combination of eigen- vectors. So start with any initial distribution ρ we like; it’s the case that n ρ= ai ρ(i) (25) i=1 8 It’s the difference between the Schr¨ odinger and Heisenberg pictures of quantum mechanics. 5 

6 .How does ρ evolve after t steps of the random walk? n ρAt = ai ρ(i) At (26) i=1 n = ai ρ(i) At (27) i=1 n = ai ρ(i) µi At−1 (28) i=1 n = ai ρ(i) µti (29) i=1 Since all of the µi < 1, except for µ1 = 1, as we step through the random walk, the distribution comes closer and closer to ρ(i) , the invariant distribution. This is ergodicity. I bring this up not only because it’s a cool bit of math, but also because it’s relevant to diffusion maps. The last few paragraphs have been all about left eigenvectors, but of course there is a close relationship between them and the right eigenvectors, which are the coordinates of the diffusion map. Recall what happens when you transpose a matrix product: (ρA)T = AT ρT (30) If ρ should happen to be a left eigenvector of A, then ρT is a right eigenvector of AT , with the same eigenvalue. Since we have arranged for our matrix A to be symmetric, A = AT , the left and right eigenvectors are the same. So we know the right eigenvalues of A — they’re the same as the left eigenvalues, starting at 1 and counting down from there towards a lowest value µn > 0. I could write the right eigenvectors as (ρ(i) )T , but that’s ugly and awkward; let’s call them f (i) instead. We’ve already seen that A1 = 1, so we’ve identified the top eigenvector: f (i) = 1. In other words, if a function is constant, then its expected value after one step is also constant. Averaging a constant just gives you back the constant. What happens to other functions, especially if we are trying to predict multiple steps ahead (i.e., averaging repeatedly)? We’ll pull the expansion-in- 6 

7 .eigenvectors trick again. n At f = At ai f (i) (31) i=1 n = ai At f (i) (32) i=1 n = ai µi At−1 f (i) (33) i=1 n = ai µti f (i) (34) i=1 → ai 1 (35) since all the eigenvalues are ≤ 1. In other words, projected far enough into the future, every function looks constant. Any irregularities in the initial function f get averaged away. This is where diffusion enters into it. In physics, diffusion refers to the process of passive drift of substances9 away from regions of high concentration and towards regions of low concentration. This is accomplished by particles of the substance taking unbiased random walks. If two adjacent regions of space start with different concentrations of the substance — say it’s higher in region 1 than in region 2 — then, even though the individual random walks are unbiased, there will be more random walkers going from region 1 to region 2 than going in the opposite direction; which will tend to reduce then difference in concentrations. More specifically, if ρ(r, t) is the density of the substance at the point r at the time t, the diffusion equation is ∂ρ ∂2ρ ∂2ρ ∂2ρ =D + 2+ 2 = D∇2 ρ (36) ∂t ∂x2 ∂y ∂z where the last equation (implicitly) defines the Laplacian ∇2 , which is some- times also written as ∆, and D is the diffusion constant.10 Remember that second derivatives are positive at minima and negative at maxima; this is say- ing that the substance, whose density is ρ, will move away from regions of high concentration and towards ones of low concentration, until it’s evenly distributed everywhere.11 On the one hand, it’s almost Biblical12 ; on the other hand, it’s 9 Or heat; which for many purposes acts roughly like a subtle fluid. 10 See the appendix. 11 Actually, precisely linear gradients are also left alone by the diffusion equation (why?). This reflects a difference between the finite-dimensional operator L for graphs, and the infinite- dimensional operator ∇2 for Euclidean space. (If this comment makes no sense, don’t worry; it’s mostly intended to reassure anyone who may have taken operator theory and started reading these notes by mistake.) 12 Isaiah 40:4, “Every valley shall be exalted, and every mountain and hill shall be laid low”. 7 

8 .the mathematical expression of what we mean when in ordinary language we say “it oozes” (Haldane, 1985, p. 33). Without going into all the ramifications of the diffusion equation (a subject that fills volumes), the important thing to notice about it is that ∇2 f says how much f changes per unit time. It takes the very simple form it does because Euclidean space is isotropic — every direction is the same as every other direction — and flat, uncurved. Now consider a diffusion process which is confined to some manifold — say some substance diffusing over the surface of a sphere. The flow can’t be given by ∇2 , since that describes inequalities of concentration in the embedding space, and we don’t care about (for instance) the fact that the concentration is > 0 on the manifold and = 0 off it — we can’t diffuse of the manifold. Rather, each manifold M has its own Laplacian operator, ∇2M , which takes into account its curvature and anisotropy.13 In fact, the connection also goes the other way: knowing a manifold’s Laplacian basically tells you its shape, because it tells you how the manifold curves. This is where we connect back, at last, to the random walk on the graph. ∇2 f says how rapidly f changes through diffusion in ordinary Euclidean space. ∇2M f says how rapidly f changes through diffusion on a manifold M. Lf says how much f changes through one step of the random walk on the graph. Suppose we build the graph by uniformly sampling points from M; as we sample more and more densely, our graph looks more and more like a discrete approximation to the continuous manifold, and so L has to encode the same geometry as ∇2M .14 So L is estimating the natural or intrinsic “shape” of the data. This remains true even if the data don’t come from uniform sampling on a manifold, but some other distribution. 2.2 Multiple Scales If A is the transition matrix for a Markov chain, then so is Am , for any integer m — it’s just taking m steps at a time, rather than one. It’s easy to show (see exercises) that Am has the same eigenvectors as A, but different eigenvalues — specifically, the eigenvalues of A, all also raised to the power m. We can therefore look at diffusion maps defined not by one step of the random walk but by m steps. The advantage of doing so is that it tends to reduce the influence of small-scale local noise. Computationally, this is because it shrinks the small eigenvalues rapidly, meaning those coordinates become less influential. 13 You might want to try to work out what the Laplacian is for a sphere, say in spherical coordinates. 14 The exact sense in which this is true is fairly subtle. On the one hand, L is an n × n matrix, so multiplying by it transforms vectors in Rn into other vectors in Rn . On the other hand, ∇2M takes functions to functions — points in RM to RM . Similarly, L represents an amount of change in one time step, while ∇2M represents a rate of change per unit time. Roughly speaking, we need to integrate ∇2M over the duration of a time-step, call it h. This 2 gives us an operator Lh defined through Lh f = eh∇M f . If we evaluate Lh f only at points on the graph, we should have Lh f ≈ Lf . See Grimmett and Stirzaker (1992) for an introduction on how discrete-time Markov chains relate to continuous-time Markov processes, and Lee and Wasserman (2008) for a precise statement of what’s going on in the present case. 8 

9 .Stochastically, it’s because, by ergodicity, where the random walk is after m steps depends less on its starting position than where it is after 1 step. Of course if m is very large, all the points blur into each other (again by ergodicity), so m is best seen as a control setting. 2.3 Choosing q It can be shown (through a very complicated argument; [[Lee and Wasserman]]) that how accurately the diffusion map reconstructs the underlying geometry depends on the ratio q i=1 µi n (37) j=1 µj where the µ are the eigenvalues of A.15 This suggests that a practical rule for choosing q is to fix a value for this ratio, and use the smallest q which achieves it. This is the default implemented in the diffuse function in the CRAN package diffusionMap. 3 What to Do with the Diffusion Map Once You Have It First, everything you might do with ordinary vector-valued data can be done with the diffusion coordinates. You can do similarity search, you can look for clusters (for instance, k-means; there’s a function to do this nicely in diffusionMap), classification, etc. We’ll be looking at regression after the midterm, and you can use the diffusion-map coordinates as the input variables in a regression. 3.1 Spectral Clustering One cute application of diffusion maps is to clustering. Remember that only the top eigenvector of A is all positive; all the other eigenvectors have both positive and negative entries. What does the difference in signs mean? To be concrete, let’s think about ρ(2) , the next-to-top eigenvector. Every point gets either a positive or a negative sign in the eigenvector. Suppose we want to start with a distribution which is concentrated solely on the points with positive sign. We can decompose any such distribution into the eigenvectors: n ρ = ρ(1) + aρ(2) + bj ρ(j) (38) j=2 but now a will be large and positive, enhancing the probability of the positive points and lowering that of the negative ones, and the bj will be small. (The bj 15 More exactly, the accuracy depends on the ratio of the sums of the population quantities of which these eigenvalues are estimates, and n in the denominator has to go to infinity. I told you it was complicated. 9 

10 .might even manage to be zero.) After one time step, this will become n (1) (2) ρA = ρ + aµ2 ρ + bj µj ρ(j) (39) j=2 . . . which looks rather like ρ: while it’s closer to being uniform than ρ was, since the eigenvalues are all < 1, it’s still concentrated on the positive points of ρ(2) , since the eigenvalues are decreasing and the bj are small. In other words, if we start with a distribution concentrated on the positive points of ρ(2) , it will tend to stay there, though diffusion will disperse it eventually. Things would work exactly the same if we concentrated the initial probability on the negative points — the only difference would be that a < 0. What about the other eigenvectors? Well, we could make a similar argument, but, because eigenvalues are smaller, initial concentrations on the positive points will diffuse away more quickly — how much more quickly will depend on the ratios of the eigenvalues. If we had to split the graph into two parts with the minimum diffusion between them, we’d split it into the points with positive and negative coordinates in ρ(2) . How does this relate to clustering? Well, suppose the data fall into two or more clusters, within which all the points are much more similar to each other than they are to outsiders. We would like to call these clusters. We could then re-arrange the kernel matrix K so it’s block-diagonal, or nearly so. Then A will also be nearly block-diagonal. Thus a random walk which starts inside one of the blocks will tend to stay inside the block for a long time. This means that all the points in a block should have the same sign in at least one of the eigenvectors. The signs of the eigenvectors, in other words, act like indicator functions, or linear combinations of indicator functions, for the clusters. In the most basic form of spectral clustering, we first divide the data into two clusters by the signs of entries in ρ(2) . Having done that, we can further sub-divide the clusters by the signs of ρ(3) and so forth. This is a top-down (divisive) hierarchical clustering. At some point it’s no longer worth it to keep splitting, and we should instead use the remaining eigenvectors as coordinates within each cluster. 3.2 Asymmetry We required the matrix K we started from to be symmetric; consequently the transition matrix A was too. In page-rank, which is otherwise similar, the transition matrix is not symmetric. How much of what we’ve done depends on symmetry, and how much would apply to any Markov chain over the data? Most of it carries through, actually. It’s no longer the case that the left and right eigenvectors are the same, but the left and right eigenvalues are. And it’s still the case that the top right eigenvector is 1, that all of the others have both positive and negative signs, etc. (Also, the non-dominant left eigenvectors all have both positive and negative signs.) It’s also the case that the signs of the entries in the non-dominant eigenvectors correspond to clusters. 10 

11 .4 The Kernel Trick Everything that went before rested on having a kernel matrix K, derived from a kernel function. Where that function came from, or what kind of data it took as arguments, was irrelevant. The data could be Euclidean vectors and the kernel the ordinary inner product, but nothing required that. So long as the kernel function was mathematically OK, nothing else mattered. This is a very powerful idea in data mining: we can split our problem into (a) finding an algorithm (for prediction or clustering or whatever) that works in terms of kernels, and (b) finding a kernel function for our representation. Then we can recycle algorithms across problems. We will return to kernel methods repeatedly, but I want to close with a small illustration of how they can be powerful, something often called “the kernel trick”.16 Suppose that we have data points x1 , x2 , . . . xn (which may be vectors or something else). Each point is represented by certain features, but we don’t think those are really the right ones for our problem. Instead, we have q different functions ψ1 , ψ2 , . . . ψq , and we really think those are the appropriate features. So we’d like to work with the vectors Ψ(x) = (ψ1 (x), ψ2 (x), . . . ψq (x)) (40) Let us say that K(xi , xj ) is the inner product of Ψ(xi ) and Ψ(xj ): q K(xi , xj ) = ψk (xi )ψk (xj ) (41) k=1 Since this is an inner product, it’s got all the properties we need for a kernel function — symmetry, forms a positive-definite matrix, etc. So if I can express my problem in terms of inner products of the transformed features Ψ, I can also write it in terms of the kernel function K on the original features. In particular, if I can find a formula for the right-hand side of Eq. 41, then I never have to calculate the functions ψ at all. In fact, I can even let q go to infinity! A concrete and very useful example is the Gaussian kernel: 1 Kh (xi , xj ) = √ exp −(xi − xj )2 /2h2 (42) 2πh2 (Sometimes you see this without the normalizing factor.) The features here are actually all the powers of xi and xj , though they’re not all equally weighted. This is in fact what people typically use with diffusion maps, rather than the Euclidean inner product, because the Gaussian distribution is preserved under diffusion. 16 For much more about kernel methods in data mining, see Shawe-Taylor and Cristianini (2004). 11 

12 .References Belkin, Mikhail and Partha Niyogi (2003). “Laplacian Eigenmaps for Dimen- sionality Reduction and Data Representation.” Neural Computation, 15: 1373–1396. Grimmett, G. R. and D. R. Stirzaker (1992). Probability and Random Processes. Oxford: Oxford University Press, 2nd edn. Haldane, J. B. S. (1985). On Being the Right Size, and Other Essays. Oxford: Oxford University Press. Edited and introduced by John Maynard Smith. Lee, Ann B. and Larry Wasserman (2008). “Spectral Connectivity Analysis.” Electronic pre-print, arxiv.org. URL http://arxiv.org/abs/0811.0121. Saul, Lawrence K. and Sam T. Roweis (2003). “Think Globally, Fit Locally: Supervised Learning of Low Dimensional Manifolds.” Journal of Machine Learning Research, 4: 119–155. URL http://jmlr.csail.mit.edu/papers/ v4/saul03a.html. Sethna, James P. (2006). Statistical Mechanics: Entropy, Order Parameters, and Complexity. Oxford: Oxford University Press. Shawe-Taylor, John and Nello Cristianini (2004). Kernel Methods for Pattern Analysis. Cambridge, England: Cambridge University Press. 12 

13 .Exercises 1. Let B be any n × n matrix. Show that if v is an eigenvector of B, then it is also an eigenvector of B2 , and of any power of B. Conclude that B and B2 have the same eigenvectors. (Hint: how many eigenvectors does each matrix have?) What happens to the eigenvalues? 13 

14 .A From Random Walks to the Diffusion Equa- tion Think of a one-dimensional space with coordinate x, and a concentration ρ(x, t) of particles of some substance at the point x at time t. Divide the axis into little intervals of width h, small enough that ρ is approximately constant over each interval. Also, let’s divide time up into intervals of duration τ , again very small. Fix a point x at time t. How does the concentration change? Assume that the particles each take independent random walks, and that the time it takes them to make one jump of the talk is much less than τ ; also that the step-size is much less than h. The expected change in position of any one particle over the time-period τ is zero, and the variance ∝ τ . Thus τD 1 τD 1 τD hρ(x, t + τ ) ≈ hρ(x, t) − hρ(x, t) + hρ(x − h, t) + hρ(x + h, t) (43) h2 2 h2 2 h2 hρ is the number of particles in the interval of length h around x (because ρ is roughly constant over the interval). The left-hand side is the new number of particles in that interval. The first term on the right hand side is the old number. The second is the number of particles which jumped out of the interval — this will shrink as τ → 0, because there’s less time for them to jump. In fact, since the variance of a random walk’s position grows ∝ t, this number really depends on τ /h2 — D is a proportionality constant. The third and forth terms count particles that jump into the interval from adjacent intervals. τD 1 τD 1 τD ρ(x, t + τ ) − ρ(x, t) ≈ 2 ρ(x, t) + 2 ρ(x − h, t) + ρ(x + h, t) h 2 h 2 h2 ρ(x, t + τ ) − ρ(x, t) D (ρ(x − h, t) − ρ(x, t)) + (ρ(x + h, t) − ρ(x, t)) ≈ τ 2 h2 ∂ρ D ∂2ρ → = ∂t 2 ∂x2 letting τ and h shrink to zero. Since D is just a proportionality constant, re-define it to absorb the factor of 1/2. This gives the one-dimensional diffusion equation, and the argument for the three-dimensional version is entirely parallel, only with more book-keeping. Of course, the idea that the particles of some substance all take independent random walks at discrete time intervals is a fable. We can pass to continuous time by making the durations and the step-sizes of the random walk smaller and smaller. The limit is a continuous stochastic process called Brownian motion (after Robert Brown, who discovered it experimentally in 1827) or the Wiener process (after Norbert Wiener, who worked out all the mathematical details in 1922). That the particles move independently of each other is a bigger assumption. It’s a reasonable approximation when they are not too highly concentrated — say molecules of perfume diffusing through the air of a room, rather than the motion of the air molecules themselves. Chapter 2 of Sethna (2006) is a good place to start. 14