Constraining Solution to be between 0 and 1.

[Dpananos]

Modelling a phenomena which only has physical interpretation between 0 and 1. The 'NonNegative` SDE option ensures the solution is positive, but additive noise my force the solution to an area in phase space which lacks interpretation.

Is there a way to constrain the solution?

[horchler]

Good question. I've worked with similar systems (N-species competitive Lotka-Volterra with additive noise).

The first thing to check is that your integration time step (specified via the tspan vector) is sufficiently small. When the noise magnitude is large or the drift term is stiff, smaller time steps may be required to not exceed limits. Ensure that the properties of your solutions remain consistent as you reduce the step size. You may see that with smaller steps, the system never exceeds the physical limits. Also make sure that you are not including a sqrt(dt) in your additive noise magnitude – this scaling is automatically applied by the solver.

In some cases, small excursions outside the limits may be acceptable. In others this may completely destabilize the system. There are a variety of ways to to solve/mitigate this issue. The 'NonNegative' option is meant as a convenience to ensure that state values can never be negative. The absolute value works similarly to a reflecting boundary condition (a max( ... , 0 ) could be used to emulate an absorbing boundary condition). It takes the Euler-Maruyama relation:

ydot(i+1) = y(i) + f(t(i), y(i))*dt + g(t(i), y(i))*dW(i)

and wraps it in a bounding function (see lines 484 and 584 of sde_euler):

ydot(i+1) = abs( y(i) + f(t(i), y(i))*dt + g(t(i), y(i))*dW(i) )

The key is that this is applied within the integration, not after it, so that the next time step won't see a state value outside of the bounds. A simple way of limiting a system between 0 and 1 would be:

ydot(i+1) = min ( max( y(i) + f(t(i), y(i))*dt + g(t(i), y(i))*dW(i), 0), 1)

To implement this in sde_euler, you could change line 484 to Yi = min(max(Yi,0),1); and line 584 to Yi(idxNonNegative) = min(max(Yi(idxNonNegative),0),1);. Note that if theses conditions are triggered often, it can alter the statistical properties of the solution. It's a good idea to examine your solution and try to adjust, the noise magnitude, parameters, and time step to avoid this in the first place if possible. I may consider adding an sdeset option for specifying custom bounding functions in the future.

There are numerous (potentially dubious) other ways to handle this situation. Many depend one the physical interpretation of the noise and what happens at the boundaries.

• Using multiplicative instead of additive noise avoids this issue around zero, but there are modeling reasons why one may prefer additive noise. More complex noise profiles that depend on the state values could be used to avoid other boundaries, provided they are acceptable to the modeler.
• One could resample from the noise until a value that doesn't exceed the limits is returned. This will of course alter the statistical properties of the solution.
• Brownian bridge interpolation could be applied to approximate solutions between time points.
• Adaptive step size SDE solvers could be used. This is an active area of research in the field (see Burrage and Burrage 2003). There are few such solvers available and they can be rather complex to implement.

Please let me know if I can further clarify anything or answer any other questions.

[Dpananos]

That is great. Thanks for the answer. I look forward to an sdeoption for this.

[alexcoca]

Hi Andreas,

I am interested in solving an SDE of the type dXt = mu(Xt,t)dt + sigma(Xt,t)dBt where mu is the gradient of the two dimensional egg-holder function and sigma(t) = sqrt(2*c/log(2+t)). However, I need to solve this given the constraints -512<=x_i<=512 (i=1,2) . I checked the sdeset and there does not seem to be an updated way of specifying constraints. I was thinking of changing the lines 484 and 584 as described above but replacing them with my constraints. Will this 'do the trick' or is a more sophisticated approach needed?

To give you a bit of background, I am looking to find the minimum of the 2-D (and 5-D) egg-holder function (number [53] here: https://arxiv.org/pdf/1308.4008.pdf) and I was thinking of comparing the simulated annealing method with solving the equation dXt = -grad(f(t)) + sqrt(2T)*dWt where Wt is a standard Brownian motion and T = c/log(2+t) where c is a constant. I am restricted to the interval (-512,512).

Many thanks,

Alex

[horchler]

@alexcoca Thanks for writing. Sorry for the delayed response.

I'm not sure I fully understand your problem. When you say -512<=x_i<=512 (i=1,2), what is x_i an how does it relate to Xt in your other equations. Is x a spatial variable - i.e., are you solving this in two dimensions in addition to time?

It would help a lot if you could formulate this as code and maybe as an ODE first so I can try to understand what you're attempting and what kind of constraints you need.

[alexcoca]

Andrew,

Thank you for your response and sorry for the rather garbled description. I'll try to make it more clear.

My task is to minimize the 5-D egg-holder function (number [53] here: https://arxiv.org/pdf/1308.4008.pdf). The minimisation is to be performed subject to the constraint that all the coordinates lie within the -512<=x_i<512 region for all i. Before I attempt the 5D problem, I decided to minimise the 2-D version that I can visualise and it's easier to work with. So yes, the answer to your question is that I am solving the equation in two dimensions to start with.

According to Geman's paper (attached), it should be possible to find a process which converges to the global minima of a function f:[a,b]^n->R by solving the SDE given by dx_t = -\grad f(x_t) dt + sqrt(2T(t)) dWt. Here x_t is a point in the domain of the function f, T(t) is a function, which is given by T = c/log(2+t) where c is a constant. The first term is just the gradient of the egg-holder function.

My question was how to capture the constrained nature of the problem, and whether simply changing 0 with -512 and 1 with 512 in the change you recommended above will do. T

The other alternative I thought of would be to include a penalty term that will increase the values of the function significantly outside the desired boundaries but is zero within the specified domain. Specifying this function is somewhat cumbersome (more so if I try to solve an SDE as opposed to using a derivative free method) so I was hoping there would be an easier way. Any suggestions/criticisms/observations/avenues to explore are welcome! I am very new to the area of stochastic optimisation so it would be a good opportunity to learn.

I shall follow up with some Matlab code implementing my idea. However, this is not at the top of my immediate priority list so it might be in the new year!

Many thanks,

Alex

Diffusions and optimization.pdf

[bykhov]

The boundaries may be easily defined by multiplication on the corresponding step functions. However, simple numerical SDE solution methods have inherent instability when discontinuity (at bounds) is introduced. Reaching the value of 512.000... is challenging. Typically, the numerical solution of equation with discontinuities requires application of implicit numerical methods (e.g. https://arxiv.org/pdf/1209.0390.pdf) rather than explicit solution generously provided by Andrew.