and r$_{j}$, j = 1 , ..., N$_{p}$ are times sampled from a Poisson process with intensity M$_{c}$ − L$_{c}$ over [0 , ˆ t ]. The variance of this estimator can be reduced by sampling N$_{p}$ from alternative distributions, as discussed in Section 4 of Fearnhead, Papaspiliopoulos and Roberts (2008). Note also that we can simulate values of ˆ u ( i ) from the proposal measure by using LEA with the modified function ˜ φ ( u ) = φ ( u ) + c ( u ), rather than use a two-step process. We now consider extending LEA to multiple dimensions. Whilst the Lam- perti transformation is usually available for univariate diffusions, it is rarely available in the multivariate case (A¨ıt-Sahalia, 2008). This is the first major restriction on the applicability of exact algorithms to multivariate diffusions. In the remainder of this article we assume that such a transformation is possible, and denote the transformed process For the proposal we use d independent Brownian motions, although we could in principle use correlated Brownian motions if required. When the do- main of the PDE of interest forms a hyperrectangle we adaptively choose the bounding interval lengths in each path segment, θ ( i ) = ( θ ( i ) $_{(1)}$, ..., θ ( i ) ( d $_{)}$), to coincide with the nearest boundary in each dimension. The methodology presented here can be adapted to any polytopal domain with appropriate rotation of the coordinate axes and selection of bounding intervals. As this algorithm shares many features with the random walk on rectangles algo- 15 d Y y s = α ( Y y $_{s}$)d s + d W W$_{s}$.