An implementation of Milstein's method for general bounded diffusions
Francisco Bernal

Slowdown of the surface diffusion during early stages of bacterial colonization
    C. Cassier-Chauvat
    December 30, 2023
===============================================================================







Despite its generality and powerful convergence properties, Milstein's method for functionals of spatially bounded stochastic differential equations is widely regarded as difficult to implement. This has likely prevented it from being utilised in applications. In this paper, we design and analyse in detail one such implementation. The presented method turns out to be on par with other, popular schemes in terms of computational cost—but with a (nearly) linear weak convergence rate under the usual smoothness requirements on coefficients and boundary. Two byproducts of theoretical interest are a new, non-standard rank-one update formula, and a connection between numerics of bounded diffusions and Eikonal equations. Three examples are worked out, confirming the accuracy and robustness of the method.
  


(Accepted for publication in the Journal of Scientific Computing, 2018)

Keywords: stopped diffusion, reflected diffusion, Feynman-Kac formula, stochastic numerics, rank-one updates to matrix decompositions, Eikonal equation. 

AMS subject classifications: 60H35, 65C30, 65C05.




§ INTRODUCTION




  
*Scope and motivation. We are interested in efficient weak schemes for Feynman-Kac functionals of stochastic differential equations (SDEs) related to general linear boundary value problems (BVPs). By efficiency, we refer to the computational cost required to bound the expected error within a given tolerance and a given confidence interval. (This depends on the weak rate of convergence w.r.t. the timestep h of the SDE scheme, and on its cost per iteration.) By general, we mean that no restrictions are posed on the coefficients and boundary beyond those required for solvability of the associated SDE (i.e. a minimum degree of regularity of both). In particular, we cover BVPs with mixed boundary conditions (BCs) (i.e. of Dirichlet type on some portion of the boundary and Robin elsewhere), for which the stochastic representation is in terms of functionals of SDEs spatially enclosed by reflecting and stopping boundaries.  

Such stochastic schemes are the basis for the pointwise solution, u( x_0,t_0), of linear BVPs in ℝ^D via the Monte Carlo method <cit.>. Besides the mathematical interest of this connection by itself, the Monte Carlo approach is computationally advantageous in many applications, v.g.   <cit.>. 
The stochastic numerics presented here can also be adapted to construct a variety of probability densities of a Wiener diffusion inside a domain enclosed by absorbing and/or reflecting surfaces, including the survival probability, the time spent close to the reflecting boundary, and first-passage times from the domain  <cit.>. 




  
*State of the art. Given a general linear BVP, it is challenging to numerically integrate its associated Feynman-Kac functionals with a satisfactory weak order of convergence w.r.t. h (δ henceforth). This is specially true when the BVP has mixed BCs—or equivalently,  when the SDE is stopped on some portion of the boundary and reflected elsewhere. The method of Constantini et al. <cit.> may be considered as the standard; it is very easy to code but has a (proven) δ of just up to one half. The difficulties pertain to the determination (in a weak sense) of the first-exit point (in the presence of stopping boundaries—see <cit.> and references therein), and of the local time (in the presence of reflecting ones). 
Other schemes have a better δ, but they are less general; the following list is not exhaustive. Gobet's half-space approximation <cit.> is  very useful in practice. He has theoretically analised the schemes for either killed (rather than stopped) or reflected diffusions; and proved that many Feynman-Kac functionals converge linearly with h—including those for homogeneous Dirichlet BCs, and those for Neumann (rather than Robin) BCs. Bossy et al. <cit.> put forward a straightforward method which has δ=1 for the Feynman-Kac functionals associated with homogeneous Neumann BCs. A relatively new approach (see <cit.> and references therein) relies on adaptive h, based on a running error estimate. It has shown very good results (including in nonconvex domains) when the diffusion is a Brownian motion, but we are not aware of extensions to more general cases. Also sticking to the Brownian motion, the method of randomisation can deal accurately with all sorts of linear BCs (see <cit.> and references therein).




Besides the method of Constantini et al., there is another one, due to Milstein, which also is completely general and has a proven O(h) weak convergence rate <cit.>: we will call it MM in short. Paradoxically (since it predates all the methods quoted above), MM has gone underreported in the literature, most likely because of a reputation of being complex to implement. (In fact, we are not aware of numerical experiments with it having been published, beyond one example in <cit.> and a few more in our review <cit.>, both on purely stopped diffusions.) For instance, MM is described as follows in <cit.>: "(...) two other weak discretization schemes (...) are considered in domains with smooth boundary (...) One of them achieves the rate of convergence h but, as pointed out by the author, is difficult to implement". Moreover, quoting from <cit.>: "An appropriate Markov chain approximation at random discretization times has been studied by Milstein (...) His procedure requires at each step near the boundary to change coordinates, and, by the way, the algorithm seems to be difficult to implement".




  
*Our contribution. The reason why MM is hard to implement is that the algorithms are incomplete. At a given point, the user is instructed to "take the ellipsoid tangent to the boundary", or to perform a given rotation or a certain matrix decomposition close to reflecting boundaries; but details as to how to are glossed over both in the seminal papers as in the later book <cit.>. In this paper, we provide a practical and theoretically sound such implementation of MM—apparently the first one.
(Henceforth, we will call it “our implementation” or "Algorithm <ref>”, to distinguish it from 
the original algorithm.[The Matlab code and data files used for this paper are available at the journal repository, and upon request from the author.])


Regarding the determination of the tangent ellipsoid (Section <ref>), we have established a novel link with an anisotropic Eikonal equation (Lemma <ref>) which allows one, theoretically, to do it exactly. Nonetheless, this approach is unpractical except with a constant diffusion matrix (as in Example <ref>). Therefore, a fast O(D^2) half-space approximation (Lemma <ref>) is also provided. For the reflections (Section <ref>), we have introduced a non-standard, non-Cholesky update formula which endowes the overall implementation with O(D^2) complexity (Algorithm <ref>). This is relevant because MM is designed for Monte Carlo simulations, which are often meant for high-dimensional problems (D≫ 1). In addition to those three main contributions, we discuss every other aspect and provide fast 
recipes for them. In sum, this paper lifts what were the main obstacles to the routinary utilization of Milstein's method for general bounded diffusions.

A preliminary, simplified version of the implementation presented here, valid only for smooth, purely stopping boundaries, was sketched in the review paper <cit.>. That algorithm, however, is now superseded by Algorithm <ref>. 

Our implementation turns out to have a cost per time step comparable (less than 2 times larger) to that of the integrator by  Constantini et al.—however, with a (nearly) linear weak convergence rate. To put this fact in perspective, the complexity of a Monte Carlo simulation (proportional to the CPU time) with error tolerance 0<ϵ≪1 drops dramatically (by a factor 1/ϵ) from δ=1/2 to δ=1 <cit.>. Substantial further gains would be possible by combining Algorithm <ref> with Multilevel <cit.>, extrapolation <cit.>, or both <cit.>. Critically, all of the previous strategies rely on the a priori knowledge of δ—thus highlighting the suitability of MM thanks to its sound theoretical foundation.  

While MM has an asymptotically proven δ=1 under the assumptions of adequate smoothness of the coefficients and of the boundary, our implementation might not exactly reproduce it in two cases:

	
  * Unless Ω has a trivial shape (or is a combination thereof), the distance map will itself be a numerical approximation—thus, possibly spoiling δ=1. (Note, however, that the same holds for any numerical scheme for bounded diffusions, not just the one presented in this paper.)

  * When the absorbing boundaries are curved and the fast recipe in Lemma <ref> for constructing the smallest tangent ellipsoid is employed (as will be typically the case). The reason is that Lemma <ref> relies on locally approximating the boundary by its tangent plane, but we have not rigourously proved convergence of δ→ 1 as h→ 0. 
  
Consequently, we claim that our algorithm has a "(nearly) linear weak convergence rate". On the other hand, we emphasize that all of the numerical experiments involving smooth boundaries which we have carried out univocally suggest δ=1, in practice.

We close the Introduction by briefly commenting on preprocessing.

The main effort which must be independently undertaken before the Monte Carlo simulation is the generation of a signed distance map for Ω. Again, we stress that this is required for any numerical scheme for bounded SDEs—although it often goes unmentioned. In the simplest cases, like a ball or a parallellepiped, an exact distance formula is available. Otherwise, we propose solving an Eikonal equation via the Fast Marching Method (check Section <ref>).

Finally, if the exact tangent ellipsoid is needed, the anisotropic Eikonal equation in Lemma <ref> must be numerically solved in advance. Example <ref>  has been crafted to illustrate a situation where this would be highly advantageous, for that problem is very hard to solve with other numerical schemes.Since preprocessing is either not particular to the algorithm presented here, or up to some point optional—and in both cases to be tackled with an independent method—its cost has not been explicitly included in the main discussion.



  
*Organisation of the paper. Section <ref> recalls the theoretical connection between linear second-order BVPs with mixed BCs, and bounded stopped/reflected SDEs. (The form of the Feynman-Kac formulas there is not the most usual one, but one tailored to MM.) In Section <ref>, the original MM is described, and the new Algorithm <ref> is listed. Implementation details are discussed around several new lemmas in Section <ref>, which is the core of the paper. Three numerical examples are worked out in Section <ref>; and conclusions are drawn in Section <ref>. To avoid clutter, all proofs of the lemmas in Section <ref> have been moved into Appendix <ref>. Finally, Appendix <ref> lists a few relevant Matlab code snippets.




§ FEYNMAN-KAC FORMULAS IN MILSTEIN'S FORM



Let D≥ 2, Ω⊂ℝ^D be a bounded domain, and Ω=Ω∪∂Ω, where the open connected set Ω is the interior of the domain and ∂Ω its boundary. Consider the linear parabolic BVP of second order with mixed BCs:


    {[ ∂ u/∂ t=  L( x,t) u + c( x,t)u + f( x,t)                        if  0<t≤ T,  x∈Ω,;                                 u= p( x)                           if  t=0,  x∈Ω,;                               u= g( x,t)                     if  0<t≤ T,  x∈∂Ω_A,;              ∂ u/∂ N= φ( x,t)u + ψ( x,t)                     if  0<t≤ T,  x∈∂Ω_R, ].


where T>0, φ( x,t)≤ 0, 
and the differential generator is given by


    L( x,t)u= 1/2∑_i,j=1^D a_ij( x,t)∂^2 u/∂ x_i∂ x_j + ∑_k=1^D b_k( x,t)∂ u/∂ x_k.


The matrix A( x,t):=[a_ij]_i,j=1^D is positive definite, and b( x,t):=(b_1,…,b_D)^T is called the drift. All of the coefficient functions in (<ref>), namely a_ij,b_i,c,f,p,g,φ and ψ are assumed continuous, and complying with the compatibility conditions at time t=0 (see <cit.> or <cit.>). The boundary is decomposed as ∂Ω=∂Ω_A∪∂Ω_R∪∂Ω_S, such that ∂Ω_A∩∂Ω_R=∂Ω_A∩∂Ω_S=∂Ω_R∩∂Ω_S=∅.
The outward[In the SDE literature, the normal is usually taken inwards. Here we follow the PDE convention.] unit normal vector N is assumed to be well defined on the boundary save perhaps on a set ∂Ω_S; 
∂Ω_A stands for the portion of the boundary (if any) where Dirichlet BCs are imposed; and on ∂Ω_R, BCs involving the normal derivative, (i.e.  N^T∇ u) hold, where ∇ u=(∂ u/∂ x_1,…,∂ u/∂ x_D)^T. (BCs involving oblique derivatives will not be considered in this paper.) Such BCs are of Neumann type iff φ=0, or of Robin type otherwise.



  
*Sufficient conditions for existence of a unique classical solution to (<ref>).
When ∂Ω_A=∅ (respectively ∂Ω_R=∅) we say the BCs are purely reflecting (resp. purely stopping), while 
when both ∂Ω_A ≠∅≠∂Ω_R, we say the BCs are mixed. If the BCs are purely reflecting (resp. purely stopping), and ∂Ω_S=∅,  theorem 2.6 (resp. theorem 2.7) in <cit.> (see also <cit.>) ensure the existence of a unique classical solution, with regularity dependent on that of the BVP coefficients and of ∂Ω. (By a classical solution, we mean that u( x,t) lives in the Hölder space C^1,2([0,T]×Ω).) In the mixed BCs case, this connection is less general and more dependent on the smoothness of the boundary <cit.>.




  
*Stochastic representation of the pointwise solution to (<ref>). Under slightly stronger conditions, the stochastic representation of the PDE with mixed BCs (<ref>) expresses its pointwise solution u( x_0,t) as the expected value of a  functional of an SDE starting at x_0 at time zero and being normally reflected on ∂Ω_R and stopped on ∂Ω_A. 

Let σ( x,t) (called the diffusion matrix) be defined by σ( x,t)σ^T( x,t)=A( x,t). (This is always possible since A is positive definite. Hence, σ≠ 0.) The following result is an extension of the well-known Feynman-Kac formulas for purely reflected (∂Ω_A=∅) and purely stopped (∂Ω_R=∅) diffusions, adapted from <cit.> and <cit.>.


	Assume that: i) a classical unique solution to (<ref>) does exist; ii) there exists a constant L>0 such that for x, y∈Ω and t∈[0,T]
	
    ||σ( x,t)-σ( y,t)||≤ L|| x- y||, 
    
    	|| b( x,t)- b( y,t)||≤ L|| x- y||;
 
	iii) ∂Ω is piecewise C^1 (i.e. C^1 save on maybe a set ∂Ω_S); and iv) either  Ω is convex, or u( x,t) can be extended to a function C^1,2([0,T]×ℝ^d). Then, for 0≤ t≤ T the pointwise solution of (<ref>) admits the following stochastic representation:
	
	
    u( x_0,t)= 𝔼[ϕ_τ]:= 𝔼[  q( X_τ)Y_τ+Z_τ ],

	
	where
	
	
    q( X_τ)={[ g( X_τ,T-τ),      if τ<T,;     p( X_T),     if τ≥ T, ].

	and the processes ( X_t,Y_t,Z_t,ξ_t) are governed by the following set of SDEs driven by a standard D-dimensional Wiener process W_t:  

    {[ d X_t= (  b( X_t,T-t) - σ( X_t,T-t)μ)dt + σ( X_t,T-t)d W_t - N( X_t)dξ_t                                                                X_0= x_0,;             dY_t= c( X_t,T-t)Y_t dt + φ( X_t,T-t) Y_t dξ_t + Y_tμ^Td W_t                                                                   Y_0=1,;            dZ_t= f( X_t,T-t)Y_t dt + ψ( X_t,T-t) Y_t dξ_t + Y_t F^Td W_t                                                                  Z_0= 0,;                                                   dξ_t=  1_{ X_t∈∂Ω_R}dt                                                                   ξ_0=0. ].

	
	Above, 1_{H} is the indicator function (1 if H is true and 0 otherwise); τ=inf_t{ X_t∈∂Ω_A} is the first exit time  (or first passage time) from Ω; which takes place at the first exit point X_τ∈∂Ω_A; and ξ_t is called the local time.

The functions μ:Ω↦ℝ^D  and F:Ω↦ℝ^D —the former a consequence of Girsanov's theorem and the latter of the expectation of Ito's integral being zero—
are at this point rather arbitrary.
However, properly choosing the function μ will later be crucial for the WoE numerical scheme.

In the remainder of the paper, F is set to zero. (See <cit.> for an application where it is not.) 

Many expected values pertaining to a population of particles can be accommodated into this framework—see <cit.>. For instance, if p=g=c=φ=ψ=0, f=1, then u( x_0,t) is the mean absorption time of a particle starting at x_0 and driven by the drift b and diffusion matrix σ in the interval [0,T]. Analogously, if p=g=c=φ=f=0, ψ=1, u( x_0,t) is the mean local time. 



  
*Elliptic equations. Equation (<ref>) can be formally 
transformed into an elliptic BVP with mixed BCs by assuming that ∂ u/∂ t=0, thus dropping the dependence on time from u and all the coefficients; letting T∞; and dropping the initial condition p. If ∂Ω_R=∅ (purely stopped diffusions / Dirichlet BCs), the stochastic representation derived from Theorem 1 still holds as long as c≤ 0 and 𝔼[τ]<∞ <cit.>. (Note that in the purely reflected case, the latter condition would be impossible.) To the best of our knowledge, there is no rigourously proved stochastic representation for elliptic BVPs with mixed BCs.  Therefore, we will assume in the remainder of this paper that, if: i) the time-independent equivalent conditions of those in Theorem 1 are in place; ii) c≤ 0; and iii) 𝔼[τ]<∞; then the time-independent equivalent representation given by Theorem 1 holds. 




§ OVERVIEW OF MILSTEIN'S METHOD AND THE PROPOSED IMPLEMENTATION (ALGORITHM 1)




 §.§ Notation


In order to approximate the SDE system (<ref>) numerically, a small, constant timestep h>0 is set, and each random realization of ( X_t)_0≤ t≤min(τ,T) is replaced by a chain (which we may call a “trajectory”) X̂_0= x_0,X̂_1,…,X̂_K, where X̂_k≈ X_kh, k=0,1,...,K, τ̂:=Kh≈τ, and t_k:=kh. The chains {Ŷ_k}_k=0^K, {Ẑ_k}_k=0^K and {ξ̂_k}_k=0^K are defined analogously. In order to lighten the notation, we drop the hats and X_k, Y_k etc. are to be understood as the discretized counterparts X̂_k, Ŷ_k, etc. unless stated otherwise. In general, functions with subindex k are evaluated at ( X_k,T-t_k), like v.g. c_k=c( X_k,T-t_k).


Vectors are by default column vectors, and ||·|| is always the Euclidean norm. For x∈ℝ^D and ∂Ω smooth, let us define: 


  * x^∂Ω=min_ y∈∂Ω|| x- y|| (i.e. the closest point on the boundary).

  * |d( x)|=|| x- x^∂Ω|| (i.e. the Euclidean distance to the boundary).

  * N( x)= N( x^∂Ω) (i.e. the normal vector pointing outwards).

  * Π( x)={ y |  N^T( x) ( y- x^∂Ω)=0} (i.e. the hyperplane tangent to ∂Ω at x^∂Ω).


We use the shorthand notation |d_k|:=|d( X_k)|, N_k:= N( X_k), and Π_k:=Π( X_k).

B( x,R) stands for a D-dimensional ball of radius R centred at x, and ∂ B( x,R) for its surface. Let M be a D× D symmetric positive definite matrix with (real) spectrum given by M m_i=μ_i m_i, μ_1≥μ_2≥…≥μ_D>0. For ρ>0, define the  D-dimensional ellipsoid centred at x and oriented according to M as

    E^ρ_M( x)= { y | ( y- x)^TM^-1( y- x)=ρ^2}.


The major semiaxis of E^ρ_M( x) has length ρμ_1^1/2 and points along m_1
[The eigenvalues and eigenvectors of M^-1 are the inverse and the same as those of M, respectively.]
, and so successively until the minor semiaxis, which runs along m_D with length ρμ_D^1/2. When M=A( X_k,T-t_k), we may simply write E^ρ_k:=E^ρ_A( X_k,T-t_k)( X_k). 

In particular, there is one value of ρ, which we call ρ, such that E^ρ_M( x) is tangent to ∂Ω. We denote it as Σ_M( x):=E^ρ_M( x)—respectively, Σ_k:=E^ρ_k. The surfaces of those ellipsoids are denoted as ∂ E^ρ_M( x), ∂ E_k^ρ, ∂Σ_M( x), and ∂Σ_k.

The eigenvalues of A( x,T-t) are λ_1( x,T-t)≥λ_2( x,T-t)≥…≥λ_D( x,T-t)>0 (or simply λ_1≥…≥λ_D if the context permits.)

The notation ω∼ W means that ω is one realization drawn from a distribution W. In particular, S_D is the uniform distribution of points over ∂ B( 0,1) (recall it is D-dimensional);  B is the distribution taking ± 1 with equal probability; and N(m,s^2) is the Gaussian with mean m and variance s^2. 





 §.§ Description of Milstein's method




  
*Remark. For the sake of clarity, let us emphasize that what we refer to as MM is the combination of two different algorithms: one for purely stopped diffusions, called “Walk on Ellipsoids” (WoE) by Milstein, <cit.><cit.>; plus the reflection off Ω_R handled according to Milstein's one-step scheme <cit.><cit.> for purely reflected diffusions. Under sufficient smoothness of the boundary and the coefficients (leading to a C^4([0,T]×Ω) solution of (<ref>)), both algorithms enjoy proven O(h) weak convergence separately. The combination of both schemes (called MM by us) is theoretically analised (and nowhere tested) in <cit.>; under the same contraints as above, it enjoys linear weak convergence rate as well. Hence, it furnishes a method for stopped/reflected SDEs and hence suitable for Monte Carlo solutions of BVPs with mixed BCs.  

Let us now explain MM, starting with its first ingredient, WoE. In contrast with Euler-like integrators,  
WoE takes bounded steps in order to avoid overshooting the boundary. 
When X_k is not interacting with it (in a sense that will be clarified in a moment), the drift is first removed from (<ref>) by setting μ=σ^-1 b. 
Then, WoE takes X_k+1= X_k+rσω, where r=√(Dh) and ω∼ S_D, whence 

    ||ω||^2 = 1 = 1/r^2(  X_k+1- X_k)^TA^-1_k(  X_k+1- X_k).


By (<ref>), the possible values of X_k+1 are distributed over ∂ E^r_k.
As long as t_k<T, WoE has the sequence { X_k}_k=1,2,... hop over ∂ E^r_0, ∂ E^r_1,… until coming close enough to the boundary (at k=i, say) that X_i iteracts with it. 
Then, one of the following happens: i) a reflection on ∂Ω_R takes place, yielding X_i+1; ii) a special ellipsoid must be taken for X_i+1; or iii) ∂Ω_A is deemed hit and the trajectory stopped.  

If X_i is so close to the absorbing boundary that E^r_i may intersect it (but still farther than r^2=Dh), WoE dictates that the tangent ellipsoid Σ_i be taken instead for X_i+1—in order to ensure that the trajectory cannot trespass the boundary. This may happen only when |d_i|<r√(λ_1( X_i,T-t_i)). Eventually, either the trajectory is stopped at the projection on ∂Ω_A (this happens when the distance to it is less than r^2), or the initial condition is hit (i.e. t_k≥ T). 

If X_i is closer than r to ∂Ω_R, the one-step scheme in <cit.> is used to handle the normal reflection on the boundary. This  involves a change of coordinates such that the new first component is pointing along - N_i.

In a nutshell, Milstein's MM is structured in the following way:

	
  * If t_k≥ T, the trajectory has 'reached' the initial condition without being stopped by ∂Ω_A. Then, read the initial condition and finish.
	
  * If not, and the trajectory is closer than r^2 to ∂Ω_A, deem it stopped. Read the Dirichlet BC and finish.
	
  * If neither of the above, and the trajectory is closer than r to  ∂Ω_R, perform the reflection, and continue.
	
  * Otherwise, hop onto the surface of E_k^r for X_k+1—if that ellipsoid is fully contained in Ω. If not, hop onto the surface of Σ_k, and continue.


As it stands, MM cannot be used, for several points must be clarified first:

	
  * the (fast) determination of the boundary data: |d( x)|, N( x), x^∂Ω, and Π( x).
	
  * the (fast) determination of λ_1( x,T-t) for general diffusions,
	
  * the (fast) determination of Σ_A( x,T-t)( x) for general domains, and
	
  * the rotation needed for handling the reflections (in a fast way).
 
By "clarified”, we mean that specific recipes must be put in place for each of the above points. Very importantly, they must be cost-efficient, since they will be performed at potentially many time steps with each of the N≫ 1 realizations involved in the Monte Carlo estimate of the expectation (<ref>). 

We prove that all of the points above can be determined at O(D^2) cost without further assumptions than those for MM. Therefore, Algorithm <ref>—which is our implementation of MM and the main result of this paper—also has an overall complexity O(D^2) per time step, which is the minimum attainable order: that of computing the matrix-vector product 
σ d W_t
[The product of a lower triangular matrix by a vector involves 1 (first row) +2 (second row) +…+D (bottom row) =D(D+1)/2 multiplications and as many additions; i.e. an O(D^2) cost.] in (<ref>). For that reason, we claim that our algorithm is "fast".
















§ DETAILED IMPLEMENTATION



The proofs of all lemmas in this Section are compiled in Appendix <ref>.



 §.§ Gershgorin estimation of the largest eigenvalue


In order to avoid overshooting the boundary with a hop from X_k, an upper bound of λ_1( X_k,T-t) is needed.
	In the event that there is no closed formula for it, accurately determining λ_1( X_k,T-t_k) (for every timestep of every trajectory in the Monte Carlo simulation) may add up to a prohibitive computational overhead. For that reason, Milstein's WoE takes a uniform upper bound 	
	λ̃_1≥max_ x∈Ωmax_0≤ t<Tλ_1( x,T-t). However, λ̃_1 itself may not be straightforward to estimate, either; or unnecessarily large, thus requiring too small an h in order to reach the required accuracy. In Algorithm <ref>, we propose substituting it by the non-uniform, sharper upper bound given below, with no detriment to δ.


	
	λ_1( x,T-t) can be bounded above at O(D^2) cost by
	
    λ_1( x,T-t)≤max_j=1,…,D∑_i=1^D|a_ij( x,T-t)|.




The cost of this approach is thus O(D) cheaper than extracting the spectrum. 




 §.§ Drift removal


In order to remove the drift from (<ref>), μ=σ^-1 b must be calculated, which may be costly, specially in high dimension—unless σ is lower triangular, where μ can be computed easily by forward substitution at cost O(D^2). 
In the solution of BVPs, σ can always be chosen lower triangular as the Cholesky matrix of A( x,t), because A is positive definite. 

On the other hand, imagine that the data is a non-lower-triangular diffusion matrix σ̂, and the point of the Monte Carlo calculation is to compute population densities such as those mentioned in Section <ref> (of the mean first exit time, etc.). Then, σ̂ should be replaced in (<ref>) by the (lower triangular) Cholesky matrix of 
σ̂σ̂^T).  In that scenario, it may be critical to obtain the lower triangular σ analytically before start, for computing the Cholesky factorization numerically at every time step involves an O(D^3) cost per step <cit.>. In the remainder of the paper, we shall assume without loss of generality that σ is lower triangular.




 §.§ Construction of the distance map



Determining the triple of boundary data {d( x), x^∂Ω, N( x)} for a point x may be time-consuming when Ω has a nontrivial shape. In general, the distance function (or distance map) inside Ω obeys the Eikonal equation in ℝ^D

    ∥∇ d( x)∥= 1,     u( x∈∂Ω)=0,      d( x∈Ω/∂Ω)<0.


With the above convention that distances are negative inside Ω, it holds

    N( x)= ∇ d( x),     x^∂Ω=  x-d( x)∇ d( x).
  

Like in all schemes for bounded SDEs, the distance map must be calculated numerically prior to the Monte Carlo simulation. As proposed in <cit.>, Sethian's Fast Marching is the method of choice.




 §.§ Construction of the tangent ellipsoid


At a given location X_k sufficiently far from the boundary, WoE draws X_k+1 from the surface of the ellipsoid E^r_k, which is inscribed in the ball B( X_k,rλ_1^1/2( X_k,T-t_k)).
When |d_k|<rλ_1^1/2( X_k,T-t_k), there is no guarantee that E^r_k does not stick out of ∂Ω.
Milstein's WoE postulates that 
the tangent ellipsoid Σ_k be taken  in that case (see <cit.>).
This prevents boundary overshoots while maximising the probability of X_k+1 being absorbed, thus leading to the least average number of hops—namely O(1/h), see <cit.>—and hence to an optimally efficient algorithm. (Even though Σ_k might be larger than E_r( X_k), the asymptotic weak convergence rate is O(h), as proven in <cit.>.)

Therefore, the determination of Σ( x) is needed for implementing Milstein's WoE. This is now formally solved by Lemma <ref> below—seemingly, a new result. 


	
	Let Ω be a closed domain in ℝ^D (not necessarily smooth).
	Then, Σ_k=E_k^Ψ( X_k)=E_k^ρ, where Ψ( x) is the solution of the anisotropic Eikonal equation 
	
    ||∇Ψ^TA_k∇Ψ||= 1 and Ψ>0  in Ω,    Ψ=0  on ∂Ω.



Lemma <ref> clarifies a fundamental issue of WoE. When A is a constant matrix, (<ref>) need be solved just once before the simulation, and then Ψ( x) will be evaluated in computing time—analogously to the Eikonal equation for the distance map. In fact, the solution to (<ref>) is required only inside a narrow shell on the inner side of ∂Ω. A Fast-Marching-like method for Lemma <ref> can solve (<ref>) only there without regard to the rest of Ω, thus cutting back on preprocessing overhead.

There are, however, two caveats to using Lemma <ref>. Numerical methods for the anisotropic Eikonal equation are
less developed (this is further commented on in Section <ref>).  
Moreover, in the event of a non-constant matrix A( x,T-t), solving (<ref>) at every X_k which needs it will in general be out of the question. 

For those reasons, we introduce the straightforward approximation of Lemma <ref>, which in Section <ref> is shown to work very well. 
The idea is to replace Σ_k by the ellipsoid tangent to the closest tangent hyperplane, which we call Σ_k.



The ellipsoid centred at X_k and tangent to the hyperplane tangent to ∂Ω 
at X_k^∂Ω 
is given by Σ_k={ y | ( y- X_k)^TA^-1_k( y- X_k)=ρ^2_k+1}, where 

    ρ_k+1=|d_k|/||σ_k^T N_k||.

When (∂Ω∩ B( X_k,|d_k|))⊂Π_k,  the ellipsoid E^r_k+1_k is fully inside Ω, where

    r_k+1=min(ρ_k+1,r).




 Lemma <ref> yields the required value r_k+1 for the next hop (check lines 18-22 in Algorithm <ref>). (Note that A_k^-1 is never used.) The cost of this construction is O(D^2), due to the product σ_k^T N_k.

When the tangency point between Π_k and Σ_k lies outside of Ω, there is a nonzero probability that X_k+1 overshoots: this is the reason of the safeguard in Algorithm <ref>, line 26. If ∂Ω is smooth at X_k^∂Ω, the portion of Σ_k sticking out of ∂Ω tends to zero as h→ 0^+ (i.e. as r→ 0^+), since Π_k→∂Ω around X_k^∂Ω. However, this may not be the case close to cusps or corners, specially if A_k is a very eccentric ellipsoid with the major semiaxis parallel to Π_k. In that worst case scenario, it is important that Σ_k shrinks with h so that the overshooting probability tends to vanish. This is the point of the rule r_k+1=min(ρ_k+1,r) in Lemma <ref>: Algorithm <ref> hops on Σ_k only if it is smaller than E^r_k (i.e. if the tangency point is nearby); otherwise it sticks to E^r_k, knowing (by Lemma <ref>) that the probability of X_k+1 overshooting goes asymptotically to zero. 
We note, however, that this is an heuristic reasoning rather than a rigourous proof that the construction given by Lemma <ref> preserves δ=1 from MM in presence of curved boundaries. 




 §.§ Change of coordinates close to the reflecting boundary
 




  
*Notation. In this subsection, column vectors in ℝ^D and ℝ^D-1 are respectively written in bold (like N_k) and with arrows (like N⃗). 
Given a D× D matrix (like σ), a starred matrix (like σ_*) denotes the submatrix obtained by removing the first column and row from the former.
 Matrix (or vector) elements superfluous for the discussion are depicted by * (like N_k^T=(*,N⃗^T)). I_S stands for the identity matrix in dimension S; and e_n is the n^th column of I_D. We drop subindex k except for X_k, X_k^∂Ω, d_k, ∇ψ_k, and N_k. All numbers are real.

When the trajectory X_k is closer than r to ∂Ω_R, coordinates are locally changed so that the new origin is X_k and the new first component points towards - N_k (i.e. inwards). The rotation is thus defined as 

    Q( X_k^∂Ω- X_k)= (-|d_k|,0,…,0)^T,

where Q is an orthogonal matrix, i.e. Q^T=Q^-1 (not unique, in general). In order to perform the rotations in an efficient way, we adapt the approach in <cit.>. 

The Givens matrix[Householder transformations could be used as well, see <cit.>.] G^θ(i,j) (where i<j) is defined element-wise as

    [G^θ(i,j)]_kl = δ_kl,  except: {[  [G^θ(i,j)]_ii= cosθ,  [G^θ(i,j)]_ij= sinθ,; [G^θ(i,j)]_ji= -sinθ, [G^θ(i,j)]_jj=  cosθ.;                       ].

Givens matrices are orthogonal and have two properties of interest to us:

	
  * Let v=(v_1,…,v_D) be an arbitrary vector. There is an angle θ( v), given by
	

    {[                   cosθ( v)=1,                    sinθ( v)=0                 if v_i=v_j=0,; cosθ( v)= v_i/√(v_i^2+v_j^2),  sinθ( v)=v_j/√(v_i^2+v_j^2),                    otherwise, ].

	
	such that the action of G^θ( v)(i,j) on a vector zeroes its j^th element, may change the i^th one, and leaves the rest unchanged:
	
    v'=G^θ( v)(i,j) v=(v_1,…,v_i-1,v'_i,v_i+1,…,v_j-1,v_j'=0,v_j+1,…,v_D)^T.

	
  * If v_i=v_j=0, and v'=G^β(i,j) v, then v_i'=v_j'=0 for any β∈ℝ.


Let v^(1) be a vector and v^(2)=G^θ_2(1,2) v^(1), where θ_2=θ( v^(1)) according to (<ref>). Then, by construction, the second element of v^(2) is zero and every other one is unchanged except for the first one. Similarly, the second and third elements of v^(3)=G^θ_3(1,3) v^(2) with θ_3=θ( v^(3)) per (<ref>) are zero. Iterating, it is clear that

    G^θ_D(1,D)⋯
    G^θ_3(1,3)
    G^θ_2(1,2) v^(1)=:G v^(1)= (β,0,…,0)^T.
 
Performing the rotations in the above order has the following properties:



Let || v^(1)||>0, and G be the matrix product of the Givens rotations G=Π_k=D^2G^θ_k(1,k) as in (<ref>), i.e. zeroing the vector elements from the second to the last, starting from v^(1). Then: (i) G_* is lower triangular; (ii) β=|| v^(1)||>0.


In order for Q N_k=(-1,0,…,0), we set v^(1)= N_k, sequentially construct G^θ_2(1,2),…,G^θ_D(1,D), and set

    Q=  -G^θ_D(1,D)⋯ G^θ_2(1,2),

and therefore

    Q^-1   = -[G^θ_2(1,2)]^-1⋯ 
    [G^θ_D(1,D)]^-1
       =
    -[G^θ_2(1,2)]^T⋯ 
    [G^θ_D(1,D)]^T
    
       =-[G^-θ_2(1,2)]⋯ 
    [G^-θ_D(1,D)].



In order to premultiply a vector, Q (or Q^T) need not be formed, but the two sets {cosθ_2,…,cosθ_D} and {sinθ_2,…,sinθ_D} are calculated and stored in advance, and later used whenever needed. To signify that the sequence of Givens rotations is performed on v, without ever constructing Q, we shall write

    Q v= -( v,θ_2θ_D),    
    Q^T v= -( v,-θ_D -θ_2).


The cost of carrying out the D-1 Givens rotations sequentially is just O(D-1), instead of O((D-1)^2) as would be the case for an explicit matrix-vector multiplication. Several vectors in Algorithm <ref> need to be rotated while computing the reflection. For instance,
(ψ̂'_1,…,ψ̂'_D)^T=Q∇ψ_k=(∇ψ_k,θ_2θ_D). 





Let us now focus on the following two parts of matrix A in the rotated frame:

    Â=QAQ^T=
    ([ Â_11=(QAQ^T)_11               *;                  Â_*= (QAQ^T)_* ]).
 
 Because A is positive definite, so are QAQ^T (since Q is orthogonal) and[If M is positive definite, M_* too, since 0< y^TM y=y⃗^TM_*y⃗, for any y=(0,y⃗^T)^T such that ||y⃗||>0.] (QAQ^T)_*, and (QAQ^T)_11>0. 
Then Â_11>0 and Â_11r^2+d_k^2>0 in Algorithm <ref>. 



  
*Updating the decomposition of the rotated submatrix. Let us now address another important computational aspect, not discussed in MM. It is required to compute the matrix-vector product Λ_*ν⃗, where  
Λ_* is defined by

    Â_*= Λ_*Λ_*^T.


Thanks to the positive-definiteness of Â_*, Λ_* could be obtained by Cholesky factorization at a cost O((D-1)^3). This is taxing if N_k is not constant on ∂Ω_R and/or in high dimension. Fortunately, the factorization A=σσ^T can be efficiently recycled into Λ_*—although, as it will be shown next, not necessarily in the standard way.
The following notation will be convenient:


    Q=
    ([   Q_11 q⃗_1^T;   q⃗_2    Q_* ]),
        σ=
    ([ σ_11 0⃗^T;   s⃗  σ_* ]).
 

Since σ is nonsingular and lower triangular by construction (see Section <ref>), σ_11≠ 0, σ_*≠ 0, and σ_* is lower triangular. On the other hand, Q_* need not be orthogonal or even nonsingular. Writing out Â= Qσσ^TQ^T gives


    Λ_*Λ_*^T= Q_*σ_*σ_*^TQ_*^T + w⃗w⃗^T,     where w⃗= σ_11q⃗_2 + Q_*s⃗.


Therefore, Λ_*Λ_*^T can be regarded as the factorization (Q_*σ_*)(Q_*σ_*)^T plus the rank-one update w⃗w⃗^T. Borrowing a standard ansatz
from <cit.>, we set


    Λ_* = Q_*σ_*(I_D-1+αz⃗z⃗^T)
  

and look for suitable α and z⃗. Inserting (<ref>) into (<ref>) and noting that (I_D-1+αz⃗z⃗^T)(I_D-1+αz⃗z⃗^T)^T= 
(I_D-1+αz⃗z⃗^T)^2=I_D-1+z⃗(α^2z⃗^Tz⃗ + 2α)z⃗^T, one has


    (α^2z⃗z⃗^T + 2α)(Q_*σ_*z⃗)(Q_*σ_*z⃗)^T= w⃗w⃗^T,


which can be readily solved by letting
α=α_+, where

    α^2z⃗^Tz⃗ + 2α=1 ⇒α_±= 
    -1±√(1+||z⃗||^2)/||z⃗||^2,


(we choose α_+ for concreteness). Then, the vector z⃗ is given by


    Q_*σ_*z⃗= w⃗.


Formulas (<ref>)-(<ref>) allow for calculation of Λ_*ν at O(D(D-1)) cost. However, Lemma <ref> below shows that this ansatz may fail. 




Let N_k^T=(N_1,N⃗^T). The update formula Q_*σ_*z⃗=σ_11q⃗_2 + Q_*s⃗ (with Q_*, q⃗_2, σ_11, σ_*, and s⃗ from (<ref>)) is inconsistent if and only if N_1=0. In that case, Q_*=0. Otherwise, there is one unique z⃗, given by

    σ_*z⃗= s⃗-σ_11/N_1N⃗.

  

We stress that Lemma <ref> is independent of the way in which the rotation is implemented. Swapping the first component for the first nozero one of N_k in (<ref>) does not help, either: the permutation matrices involved induce structural changes (the analogous of σ_* is no longer lower triangular, for instance), with the result that Milstein's formulas in Algorithm <ref> would have to be reworked. 

When the first component of N_k is zero, Λ_* can still be computed at O((D-1)^2) cost by  using Algorithm  in LINPACK <cit.> (implemented as  in Matlab).
This exploits Lemma <ref>: since Q_*σ_* is lower triangular, Λ_* can be seen as a rank-one update of a Cholesky factorization. Despite the fact that Q_*σ_*σ_*^TQ_*^T is only positive semidefinite if N_1=0,  would also work in that case—thanks to the fact that one Cholesky matrix, namely Q_*σ_*, is available in the first place. However, using  would involve forming the full matrices. 

Instead, we put forward the following analytical approach, specifically tailored to the case N_1=0. It relies on the following decomposition:


    Λ_*Λ_*^T=
    
    [Q_*σ_* | w⃗]
    [[ σ_*^TQ_*^T;           ;       w⃗^T ]]=
    [Q_*σ_* | w⃗]
    ZZ^T
    [[ σ_*^TQ_*^T;           ;       w⃗^T ]],


where Z is an orthogonal matrix of order D. If it can be chosen such that


    Z^T
    [[ σ_*^TQ_*^T;           ;       w⃗^T ]]=
    [[    H;     ; 0⃗^T ]],


then Λ_*=H^T. It turns out that Z can be found analytically, leading to a closed formula for Λ_*.



Let N_k^T=(0,N⃗^T); Q_*, q⃗_2, σ_11, σ_*, and s⃗ from (<ref>),
and w⃗ from (<ref>). Then,

    Z=
    [[ I_D-1-1/||r⃗||^2r⃗r⃗^T              r⃗/||r⃗||;                       ;            r⃗^T/||r⃗||                      0 ]],

where σ_*r⃗=N⃗. Furthermore,  Λ_*=Q_*σ_*+1/||r⃗||w⃗r⃗^T.




  
*Remark.  Λ_* in Lemma <ref> is not necessarily triangular—but note that this was not required, anyways (check <cit.>, where Λ_* is called λ).


Let v⃗ be an arbitrary column vector in ℝ^D-1. Observe that

    Q([  0; v⃗ ])=
    ([     *; Q_*v⃗ ])=    -( 
    ([  0; v⃗ ])
    ,θ_2→θ_D ),
    
    Q^T([  0; v⃗ ])=
    ([       *; Q_*^Tv⃗ ])=    -( 
    ([  0; v⃗ ])
    ,-θ_D→ -θ_2 ).


Let [QM]_n be the n^th column of QM, where M is a square matrix of order n. Then, [QM]_n can be calculated as [QM]_n=(QM) e_n=-(M e_n,θ_2→θ_D). Furthermore, the first column of  matrix Â is QAQ^T e_1=-Q(σσ^T N_k), i.e.

    ([ Â_11;    ⋮; Â_D1 ])=
    Â e_1
    = 
    (
    σσ^T N_k
    ,θ_2θ_D).


The concrete calculation of Λ_*ν with cost O(D(D-1)) is listed as Algorithm <ref>. Complexity is dominated by the product σ_*ν⃗ and the forward substitutions.









§ NUMERICAL EXPERIMENTS


In this section, we report numerical results obtained with Algorithm <ref>—henceforth, they are labelled as "MM”. The algorithm of Constantini et al. has also been used for comparison (“REF”). (Specifically, the version with pseudonormal variables, which has a weak convergence rate O(h^1/2-ϵ),ϵ>0 <cit.>.) Both codes have been written in Matlab, and run on a laptop. 
In that case (even though Algorithm <ref> is written in sequential form), the codes should be fully vectorized, which is critical for speed (check appendix <ref>). 

The Matlab code and data files used for this paper are available at the journal repository, and upon request from the author. 

Let
ϕ_h^(1),…,ϕ_h^(N) be N iid numerical approximations to ϕ in  (<ref>). The numerical approximation to the Feynman-Kac functional is 

    u_h,N( x_0,t)= 1/N∑_j=1^N ϕ_h^(j).

Asymptotically (i.e. as N→∞ and h→ 0^+),  ε_N,h∼ N(Ch^δ,𝕍[ϕ_h]/N), where ε_N,h is the root mean-square (RMS) error of ϕ_h, C>0 is a constant and δ>0 is the weak convergence rate of the scheme <cit.>. In all the ensuing experiments, the pointwise exact solution u_ex( x_0,T)=𝔼[ϕ]≠ 0 is known. Then, after setting an accuracy goal ε, N is chosen accordingly, i.e. such that 2√(𝕍[ϕ_h]/ N)= .20×ε. The idea is that asymptotically the RMS error carries less than a 20% statistical error with a large probability. 
However, for better assessment and comparison across experiments, the relative error ϵ_h^rel:=|1-u_h,N( x_0,T)/u_ex( x_0,T)| is reported instead. 



  
Remark. In practice, one starts simulating trajectories j=1,2,..., replaces 𝕍[ϕ_h] by the sample estimate so far, V_j, and stops as soon as √(V_j/j)≤ε/10.




 §.§  Example I (three-dimensional)


We take this from <cit.>. The coefficients of (<ref>) are:
c=0, b=(y,z,x)^T, φ=-|| x||^2, 

    σ(x,y,z)=([       √(1+|z|)              0              0;    1/2√(1+|x|) √(3/4)√(1+|x|)              0;              0    1/2√(1+|y|) √(3/4)√(1+|y|) ]).


By choosing as exact solution u_ex( x,t)=xyz (independent of time), and g=p=u_ex, the remaining coefficients ψ and f are derived from (<ref>), i.e. ψ= N^T∇ u_ex-φ u_ex, and 
f= y^2z + z^2x + x^2y + 1/2√(1+| z |)√(1+| x |)+√(3)/2x√(1+| x |)√(1+| y |).

As in <cit.>, we set Ω=B( 0,1)  and x_0=(.56,.52,.30)^T; thus 
u_ex( x_0,T)≈.08736. We consider three different sets of BCs: purely absorbing, purely reflecting, and mixed. In the latter case, the hemisphere with z<0 is absorbing and the other one, reflecting. (Note that the reflection is not conormal.) The convergence of 
ϵ_h^rel w.r.t. h is shown in Table <ref>, along with the results with REF. Regardless of the BCs, MM is both more accurate and has a faster convergence rate (estimated by least-squares regression), in fact very close to the theoretical value δ=1. 





In many applications boundaries are not so smooth as on a sphere; we also solved this problem in the box Ω=[-√(2)/2,√(2)/2]^3. We only show results for the purely absorbing and purely reflecting cases; see Table <ref>. Computational times have been included.

Since the boundary is nonsmooth now (due to the corners), accuracy and convergence rate are bound to deteriorate. In particular, theoretical rates no longer apply. Nonetheless, MM is still the more efficient integrator. For instance—acording to Table <ref>—REF took 340 s. to attain a relative error of .12722, while MM took just 65.7 s. for .12573 (with the same confidence interval). We stress that this is the case where MM performs worst (purely reflecting nonsmooth boundary).






 §.§ Example II (arbitrary-dimensional)


This problem features crossed second derivatives, oscillating coefficients and solution, and can be posed in any dimension. The diffusion matrix is

    σ=([ 1 0 … 0; 1 1 … 0; ⋮ ⋮ ⋱ ⋮; 1 1 … 1 ])
    ⇒σσ^T= A= 
    ([ 1 1 1 … 1; 1 2 2 … 2; 1 2 3 … 3; ⋮ ⋮ ⋮ ⋱ ⋮; 1 2 3 … D ]).

The exact solution is u_ex=cos∑_i=1^D x_i. The coefficients are: b=(sinπ x_1,…,sinπ x_D)^T, c=0, φ=-1, and g=u_uex (f and ψ can be derived from (<ref>)). We take T=∞ (an elliptic PDE, so that no initial condition is needed), and Ω=B( 0,1) in ℝ^D.





In Table <ref>, we consider the same mixed BCs as before (i.e. reflecting on the upper hemisphere). The solution with MM at x_0=(-0.7854,0,…,0)^T is evaluated for increasing D. For comparison, the solution with REF in D=5 is also included. Times are not included, but let us give an example: REF took 3119 s. for a relative accuracy of .01248, while MM (always in D=5) took 50.9 s. for ϵ_h^rel=.01320. As usual, the difference grows with increasing accuracy (as O(1/ϵ_h^rel), see <cit.>) due to the better δ of MM. 

We use Example II to compare the cost per iteration of MM and REF. Specifically, we solve it in B(0,1) with T=1, same x_0, N=10^5 (fixed),  h=.001 and purely reflecting BCs (so that the number of steps is roughly constant). The total times taken by MM and REF are shown on the left of  Figure <ref>, for D=3,…,8. This indicates that the mean cost per iteration of MM is less than twice as with REF.

All numerical experiments so far have made use of Gershgorin's circles theorem (Section <ref>) for estimating λ_1. We checked that doing this does not have a discernible impact on accuracy (compared with using the exact λ_1), while it cuts back on computational overhead (by about 15%-20%, in these examples).





 §.§ Example III (two-dimensional with sensitive corners)
 
The purpose of this problem is to investigate the connection between MM and the anisotropic Eikonal equation established in this paper. We pick T=∞,  x_0=(.823,-.875)^T and purely absorbing BCs in Ω=[-1,1]^2. The coefficients are:

    A=
    ([     8 -2.71; -2.71     1;       ])
    =σσ^T
        ⇒    σ≈([  2.8284       0; -0.9581  0.2863;         ]),

c=0, b=(sinπ x,sinπ y)^T, g=u_ex, and f is derived from the exact solution,

    u_ex(x,y)=1/x^2+y^2-R.

The point of (<ref>) is that, as R→ 2, the exact solution tends to infinity on the four corners, substantially compounding the difficulty of integrating the corresponding SDEs in the domain with nonsmooth boundary. 
Problem III is solved with MM and REF, as before, plus a new version of MM where, instead of using Lemma <ref> for Σ_k, Σ_k is calculated exactly according to Lemma <ref>—we call this integrator “MM+”. In order to do so, the anisotropic Eikonal equation  

    ||∇Ψ^TA∇Ψ||= 1 and Ψ>0  in Ω,    Ψ=0  on ∂Ω,

is solved before the Monte Carlo simulation. (Importantly, note that thanks to the fact that A is constant, only one such equation must be solved.) We do so using the Fast Sweeping scheme described in <cit.>, and store the numerical solution in a 251× 251 grid, from where Ψ( X_k) is later interpolated “on the fly”.  (Lines 19-21 in Algorithm <ref> must then be replaced by r_k+1=Ψ( X_k).)

The contours of Ψ are sketched on Figure <ref> (right). Figure <ref> shows the convergence of ϵ_h^rel for increasing difficulty: R={2.1,2.01,2.0001}. Note that MM+ converges faster than MM, and both versions can tackle this problem much better than REF—which in the hardest case R=2.0001,  essentially cannot cope.



§ CONCLUSIONS


Numerical experiments 
show a very satisfactory performance of Algorithm <ref>, in terms of accuracy, speed, and weak convergence rate.  
We stress that it has a favourable complexity, rendering it suitable for high-dimensional problems.

The connection, established in this paper, with the anisotropic Eikonal equation, paves the way for the efficient numerical integration of general absorbed diffusions in domains with corners. (The only other method with this capability known to the author is the Walk on Rectangles algorithm <cit.>, which is nonetheless restricted to the Brownian motion and some variations thereof with constant coefficients.) On the other hand, solving the anisotropic Eikonal equation is nontrivial and the object of intense current research. Tailoring existing solvers for the anisotropic Eikonal equation to the particular features of the algorithm presented here is left as future work.

Finally, we point out that the most promising direction for further development of this algorithm is producing a Multilevel formulation, along the lines of <cit.>.

  




§ ACKNOWLEDGEMENTS

The author thanks J. A. Acebrón for proposing this work and  helpful discussions throughout. Portuguese FCT funding is gratefully acknowledged. 



§ PROOFS OF LEMMAS IN SECTION 4


Proof of Lemma <ref>. This is a direct application of Gershgorin's circles theorem: for any eigenvalue λ of any D× D matrix M, there is 1≤ i ≤ D such that

    |λ-M_ii|≤∑_i≠ j|M_ij|.

Since A is positive definite, λ>0. By the triangle inequality, |λ-a_ii|≥ |λ|-|-a_ii|=λ-|a_ii|, and hence

    λ_1( x,T-t)≤max_i=1,…,D∑_j=1^D|a_ij( x,T-t)|.

The cost is that of adding up D elements for each of the D rows, and comparing those D values to find the largest one. 


Proof of Lemma <ref>. Let M be a D× D symmetric positive definite matrix, with spectrum M m_i=μ_i m_i,  where μ_1≥…≥μ_D>0 and || m_i||=1. The (positive) viscosity solution of the anisotropic Eikonal equation

    ||∇Ψ^TM∇Ψ||= 1  in Ω,    Ψ=0  on ∂Ω,

at x∈Ω, namely Ψ( x), can be interpreted as the arrival time at x of a monotonic front which marches anisotropically with velocity μ_i along the direction m_i (see <cit.> and references therein). 

By symmetry, it is also the arrival time on ∂Ω of a front marching outwards from x with the same velocities in the opposite directions (i.e. μ_i along - m_i).

Then, at time t'>0, the front stemming from x is the ellipsoid with semiaxes of length μ_1t',…,μ_Dt' along the orthogonal directions m_1,…, m_D. Since M and M^-1 have the same eigenvectors
and inverse eigenvalues, that ellipsoid is the locus of ( y- x)^TM^-1( y- x)=t'^2. 

In particular, at time t”=Ψ( x), the ellipsoidal front arrives on the boundary for the first time, meaning that it is tangent to it. Therefore,

    Σ_M( x)= { y | ( y- x)^TM^-1( y- x)=Ψ^2( x)}.
 
Setting M=A_k and Σ_k=E^|Ψ( X_k)|_A( X_k,T-t_k)( X_k) yields the desired tangent ellipsoid.

Proof of Lemma <ref>. Without loss of generality, let us take the origin at X_k. Set m( y,ρ)=  y^TA^-1_k y-ρ^2.  
The sought-for ellipsoid is then Σ_k={ y | m( y,ρ_k+1)=0}. The tangency point y_0 belongs both to Σ_k and to Π_k, so that for some q≠ 0,

    ∇ m( y_0,ρ_k+1)=2A^-1_k y_0=q  N_k.

Since A_k≠ 0, (<ref>) implies

    y_0=q/2A_k N_k⇒{[ q^2/4 N_k^TA_k^TA^-1_kA_k N_k= ρ^2_k+1,
     and;                  N_k^T y_0=q/2 N_k^TA_k N_k. ].

Moreover, A_k=A^T_k and hence

    ρ_k+1=  N_k^T y_0/√( N_k^TA^T_k N_k)=
    |d_k|/||σ^T N_k||.


The condition (∂Ω∩ B( X_k,|d_k|))⊂Π_k is equivalent to ∂Ω being planar inside B( X_k,|d_k|). Let us assume it is. Then, E^r_k sticks out of ∂Ω iff it sticks out of Π_k. If it doesn't, then ρ_k+1≥ r. If it does, then Σ_k=E^ρ_k+1_k is concentric to and inside of E^r_k, so that ρ_k+1<r. In both cases, E^r_k+1( X_k) given by rule (<ref>) is inside Ω. 

Proof of Lemma <ref>. In preparation, let θ∈ℝ, 1<j,k≤ D, and w=G^θ(1,k) e_j . Then, 

    if k≠ j and j>1⇒ G^θ(1,k) e_j=  e_j.

To see this, simply note that w_1=( e_j)_1cosθ+( e_j)_ksinθ and w_k=-( e_j)_1sinθ+( e_j)_kcosθ are both zero if j>1 and j≠ k, while the other elements of w are unaffected by the rotation and thus are the same as in e_j.

Assume now 1<i<j≤ D. Then,

    G_ij   =  e_i^TG e_j 
       =  e_i^T G^θ_D(1,D)⋯ G^θ_i+1(1,i+1) G^θ_i(1,i) ⋯ G^θ_2(1,2)  e_j 
       =([G^θ_i+1(1,i+1)]^T⋯[G^θ_D(1,D)]^T e_i)^T (G^θ_i(1,i) ⋯ G^θ_2(1,2)  e_j)
       =(G^-θ_i+1(1,i+1)⋯ G^-θ_D(1,D) e_i)^T e_j
       = e_i^T e_j=0.


In (<ref>), we have used the facts that: [G^θ(m,n)]^T=G^-θ(m,n); that Π_k=i^2G^θ_k(1,k) e_j= e_j because the interval 2,…,i does not include j, by assumption; and analogously i+1,…,D does not include i. 
Since G_1<i≤ D,1<j<i=0, G_* is lower triangular, thus proving (i).

On the other hand, the upper row of G is not zero, in general. In fact, working out those entries (which is tedious but straightforward), one has:

    G_11   = cosθ_2cosθ_3cosθ_4⋯cosθ_D,
    
    G_1k   =sinθ_kcosθ_k+1cosθ_k+2⋯cosθ_D, (2≤ k≤ D).




Let us now prepare for (ii). First, let v be a vector and v'=G^θ( v)(1,k) v. By (<ref>),

    v'_1=√((v_1^(1))^2+(v_k^(k))^2)≥ 0.


Consequently, in the sequence v^(k)=G^θ_k(1,k) v^(k-1) with k=2,…,D, it holds that v_1^(2)≥ 0,…,v_1^(D)≥ 0. By (<ref>),  cosθ_3,…,cosθ_D are thus all nonnegative. (Only cosθ_2<0 iff v_1^(1)<0). 

Clearly, β^2=( v^(1))^TG^TG v^(1)=
( v^(1))^T v^(1)= || v^(1)||^2>0, by hypothesis. 
Also, it is always possible to pick n such that v^(1)_n≠ 0. Then, v_n^(1)= e_n^T v^(1)=
 e_n^TG^TG v^(1)=
 e_n^TG^T(β,0,…,0)^T=
β G_1n. This prevents G_1n from being zero, so that β=v_n/G_1n. Moreover, by the previous discussion and (<ref>), G_1n≠ 0 implies that cosθ_n+1⋯cosθ_D>0. 

Therefore, it is clear that, if n=1, sign(G_11)=sign(cosθ_2)=sign(v_1^(1)), by (<ref>) and (<ref>).


If 1<n≤ D, sign(G_1n)=sign(sinθ_n)=sign(v_n^(n-1)), by (<ref>). Let us show that v_n^(n-1)=v_n^(1). This occurs because v^(k)=G^θ_k(1,k) v^(k-1) preserves all but the first and k^th elements of v^(k-1), so that if n>k, then v_n^(k)=v_n^(k-1). Iterating, it is clear that v^(n-1)=G^θ_n-1(1,n-1)⋯ G^θ_2(1,2) v^(1) has preserved the elements v^(n-1)_n,…,v^(n-1)_D from v^(1). Then, by (<ref>),  sign(sinθ_n)=sign(v_n^(n-1))=sign(v_n^(1)).

Summing up, we have proved that sign(β)= sign(v^(1)_n/G_1n)>0. In fact, β=+|| v^(1)||, which is (ii). 



Proof of Lemma <ref>.  We recall that σ_*≠ 0 and σ_11≠ 0 because σσ^T is positive definite. By definition Q(N_1,N⃗^T)^T=(-1,0⃗^T)^T and hence

    Q_11N_1+q⃗_1^TN⃗=-1,


    N_1q⃗_2+Q_*N⃗=0⃗.

Writing out QQ^T=I_D and Q^TQ=I_D in terms of the blocks defined in  (<ref>) yields

    Q_11^2 + q⃗_1^Tq⃗_1= 1,


    Q_11^2 + q⃗_2^Tq⃗_2= 1,


    -Q_11q⃗_2 = Q_*q⃗_1,


    -Q_11q⃗_1 = Q_*^Tq⃗_2,


    q⃗_2q⃗_2^T + Q_*Q_*^T= I_D-1,


    q⃗_1q⃗_1^T + Q_*^TQ_*= I_D-1.





	
	
  * Let us start by the 'if' part. Assume that N_1=0. 
	
	Since || N_k||=1, then N⃗_k≠0⃗. From (<ref>), Q_*N⃗=0⃗, and therefore Q_*=0. Applying Sylvester's determinant identity to (<ref>) gives:
	
    (Q_*Q_*^T)= (I_D-1-q⃗_2q⃗_2^T)= 1-q⃗_2^Tq⃗_2  ⇒ q⃗_2^Tq⃗_2=1.

	
	Then, by (<ref>), Q_11=0. By (<ref>), this means that q⃗_2^TQ_*=0⃗^T, so that q⃗_2 is not in the range of Q_* (i.e. it is orthogonal to the column space of Q_*^T). On the other hand, Q_*(σ_*z⃗) has a projection on q⃗_2, since q⃗_2^TQ_*(σ_*z⃗)=σ_11≠ 0. Consequently, there cannot be such σ_*z⃗ if N_1=0.
	
	
  * Assume now that the system is inconsistent. Since σ_11≠ 0 and Q_*s⃗ in w⃗ clearly belongs to the range of Q_*, it follows that q⃗_2≠ 0 does not. This implies q⃗_2^TQ_*=0⃗^T, and by (<ref>), N_1=0.  This proves the “only if” part.
	
	
  * Finally, assume there is a solution. It has been proved that then N_1≠ 0 and Q_11≠ 0. By Sylvester, (Q_*Q_*^T)=1-q⃗_2^Tq⃗_2=Q_11^2, whence Q_*≠ 0.  This means that σ_*z⃗, and thus z⃗, is unique.
	
	If N_1≠ 0, then -σ_11/N_1Q_N⃗=σ_11q⃗_2 by (<ref>). It follows that Q_*( s⃗-σ_11/N_1N⃗)=w⃗. 
	





Proof of Lemma <ref>. Let us try to fit (<ref>) with


    [[    Φ   r⃗; r⃗^T    0 ]]
    [[ σ_* ^TQ_*^T;            ;        w⃗^T ]]
    =
    [[ Φσ_*^T Q_*^T + r⃗w⃗^T;       r⃗^Tσ_*^T Q_*^T ]].


In order for the bottom row of the rightmost matrix above to vanish, Q_*σ_*r⃗=0⃗.  Since σ_*≠ 0, this has a nontrivial solution iff Q_*=0, i.e. iff N_1=0, by Lemma <ref>. In that case, by (<ref>), σ_*r⃗=λN⃗, with λ∈ℝ. (For definiteness, we pick λ=1.)

The second requirement for a candidate Z is that it be orthogonal. By symmetry of the righmost matrix in (<ref>), it suffices to check that Z^TZ=I:


    Z^TZ=
    [[ ΦΦ^T+r⃗r⃗^T         Φr⃗;     r⃗^TΦ^T      r⃗^Tr⃗ ]]
    =
    [[ I_D-1    0⃗;  0⃗^T     1 ]]
    ⇔{[ ΦΦ^T=I_D-1-r⃗r⃗^T               (i);            Φr⃗=0⃗              (ii);          ||r⃗||=1             (iii) ].


Assume for the time being that ||r⃗||=1. Let us take Φ=I_D-1-r⃗r⃗^T. Then
Φ=Φ^T and ΦΦ^T=Φ^2=I_D-1-2r⃗r⃗^T+r⃗(r⃗^Tr⃗)r⃗^T=I_D-1-r⃗r⃗^T=Φ, meeting (i). For (ii), note that Φr⃗=(I_D-1-r⃗r⃗^T)r⃗=r⃗-(r⃗^Tr⃗)r⃗=0⃗. Finally, condition (iii) is met simply by substituting r⃗/||r⃗|| for r⃗ (this makes the result independent of λ.)


Therefore, Z in (<ref>) is orthogonal, and yields Λ_*=Q_*σ_*(I_D-1-1/||r⃗||^2r⃗r⃗^T)+1/||r⃗||w⃗r⃗^T. Replacing r⃗ by σ_*^-1N⃗ and using the fact that Q_*N⃗=0⃗ finishes the proof. 




§ MATLAB CODE SNIPPETS


Vectorized functions are given below for  (Section <ref>), forward subtitution with a lower triangular system, and distribution B (Algorithm <ref>). 













10




Bayer2010
C. Bayer, A. Szepessy and R. Tempone, Adaptive weak approximation of reflected and stopped diffusions. Monte Carlo Methods App. 16, 1-67 (2010).

IterPDD
F. Bernal and J.A. Acebrón, A multigrid-like algorithm for probabilistic domain decomposition. Comput.  Math. App. 72(7) 1790–1810 (2016).


OurReview
F. Bernal and J.A. Acebrón, A comparison of higher-order weak numerical
schemes for stochastic differential equations in bounded domains.
Comm. Comput. Phys. 20(3) 703–732 (2016).

Bernal2014
F. Bernal, J.A. Acebrón and I. Anjam, A Stochastic Algorithm Based on Fast Marching for Automatic Capacitance Extraction in Non-Manhattan Geometries.
SIAM J. Imaging Sci. 7(4), 2657-2674 (2014).



BossyGobetTalay
M. Bossy, E. Gobet, and D. Talay, A symmetrized Euler scheme for an efficient approximation of reflected diffusions, J. App. Probability 41(3) 877–889 (2004).

Buchmann_Petersen_SIAM
F.M. Buchmann and W.P. Petersen, An Exit Probability Approach to Solving High Dimensional Dirichlet Problems. SIAM J. Sci. Comput. 28, 1153-1166 (2006). 


Cao2007 
N. Cao, M. Ortner, A. Nehorai,  Solutions for diffuse optical tomography using the Feynman-Kac formula and interacting particle method, Proc. SPIE BIOS'07, San Jose, CA (USA) Vol. 6434 (2007) 643402-64313.

Constantini1998
C. Constantini, B. Pacchiarotti and F. Sartoretto, Numerical approximation for functionals of reflecting diffusion processes. SIAM J. Appl. Math. 58, 73-102 (1998).

Deaconu2006
M. Deaconu and A. Lejay, A Random Walk on Rectangles Algorithm, Methodol. Comput. Appl. Probab. 8, 135-151 (2006).

Dongarra_LINPACK
J. J. Dongarra,  J. R. Bunch, C. B. Moler, and G. W. Stewart. LINPACK Users' Guide. SIAM, Philadelphia (1979). 

Freidlin_Book
M. Freidlin, Functional Integration and Partial Differential Equations. Annals of Mathematics Studies, vol. 109, Princeton University Press (1985).


Giles_y_yo
M.B. Giles and F. Bernal, Multilevel simulations of expected exit times and other functionals of stopped diffusions. SIAM / ASA J. Uncertainty Quantification, 6(4), 1454–1474 (2018). 


Gill_matrix_updates
P.E. Gill, G.H. Golub, W. Murray, and M. Saunders, Methods for modifying matrix factorizations.
Math. Comput., 28 505-535, (1974).


Gobet   Menozzi_2010
E. Gobet and S. Menozzi, Stopped diffusion processes: overshoots and boundary correction. Stoch. Proc. Appl. 120, 130-162, (2010).


Gobet2001_BB
E. Gobet, Euler schemes and half-space approximation for the simulation of diffusions in a domain. ESAIM: Prob. Stat. 5, 261-297 (2001).

Gobet2016 
E. Gobet, Monte Carlo Methods and Stochastic Processes. Chapman & Hall/CRC (2016).

Rusos 
O.A. Ladyženkaja, V.A. Solonnikov, N. Ural'ceva, Linear and Quasi-linear
Equations of Parabolic Type. Trans. Math. Monographs 23, AMS, Providence, RI (1968).



Pages
V. Lemaire and G. Pagès, Multilevel Richardson-Romberg extrapolation. Bernoulli
23(4A), 2643–2692 (2017).

Lepingle1995 
D. Lépingle, Euler scheme for reflected stochastic differential equations, Math. Comput. Simul. 38 119-126 (1995).

Licata   Grill 
N.A. Licata, S.W. Grill, The first passage problem for diffusion through a cylindrical pore with sticky walls, Euro. Phys. J. E 30(4)  439–447 (2009).


Lieberman1986 
G.M. Lieberman, Mixed boundary value problems for elliptic and parabolic differential equations of second order,
J. Math. Anal. App. 113(2), 422–440 (1986).


Maire   Simon
S. Mairé and M. Simon, A partially reflecting random walk on spheres algorithm for electrical impedance tomography, J. Comput. Phys., 303 413–430 (2015).


Sara
S. Mancini, F. Bernal and J.A. Acebrón, An efficient algorithm for accelerating Monte Carlo approximations of the solution to boundary value
problems. J. Sci. Comput., 66(2), 577–597 (2016).

Mascagni2002 
M. Mascagni, N.A. Simonov, Monte Carlo methods for calculating
some physical properties of large molecules, SIAM J. Sci. Comput. 26, 339-357  (2004).


Milstein_Neumann
G.N. Milstein, Application of the numerical integration of stochastic equations for the solution of boundary value problems with Neumann boundary conditions, Theor. Prob. Appl. 41, 170-177 (1996).

Milstein_WoE 
G.N. Milstein, Weak Approximation of a Diffusion Process in a Bounded Domain. Stoch. Stoch. Rep., 62 147-200 (1997).

Milstein_Tretyakov_Book
G.N. Milstein and M.V. Tretyakov, Stochastic Numerics for Mathematical Physics. Springer, Berlin (2004). 



Miranda1970 
C. Miranda, Partial Differential Equations of Elliptic Type, 2nd ed. Springer-Verlag, Berlin/Heidelberg/New York (1970).


Redner2001
S. Redner, A Guide to First-Passage Processes. Cambridge University Press (2001).

Schwabedal_PRL
J.T.C. Schwabedal and A. Pikovsky, Phase Description of Stochastic Oscillations, Phys. Rev. Lett. 110(20), 204102 (2013).

Schwenke2012
M. Schwenke, A. Hennemuth, B. Fischer, and O. Friman, A Novel Anisotropic Fast Marching Method and Its Application to Blood Flow Computation in Phase-Contrast MRI. Methods Inf. Med. 51(5) 423–428 (2012).


Fast_Sweeping
Y.R. Tsai, L.T Cheng, S. Osher, and H.K. Zhao, Fast sweeping algorithms for a class of
Hamilton-Jacobi equations. SIAM J. Numer. Anal. 41(2), 673–694 (2003).




