Multilevel and multi-index sampling methods with applications - Lecture 1: Adaptive strategies for Multilevel Monte Carlo
Tempone, Raul
We will first recall, for a general audience, the use of Monte Carlo and Multi-level Monte Carlo methods in the context of Uncertainty Quantification. Then we will discuss the recently developed Adaptive Multilevel Monte Carlo (MLMC) Methods for (i) It Stochastic Differential Equations, (ii) Stochastic Reaction Networks modeled by Pure Jump Markov Processes and (iii) Partial Differential Equations with random inputs. In this context, the notion of adaptivity includes several aspects such as mesh refinements based on either a priori or a posteriori error estimates, the local choice of different time stepping methods and the selection of the total number of levels and the number of samples at different levels. Our Adaptive MLMC estimator uses a hierarchy of adaptively refined, non-uniform time discretizations, and, as such, it may be considered a generalization of the uniform discretization MLMC method introduced independently by M. Giles and S. Heinrich. In particular, we show that our adaptive MLMC algorithms are asymptotically accurate and have the correct complexity with an improved control of the multiplicative constant factor in the asymptotic analysis. In this context, we developed novel techniques for estimation of parameters needed in our MLMC algorithms, such as the variance of the difference between consecutive approximations. These techniques take particular care of the deepest levels, where for efficiency reasons only few realizations are available to produce essential estimates. Moreover, we show the asymptotic normality of the statistical error in the MLMC estimator, justifying in this way our error estimate that allows prescribing both the required accuracy and confidence level in the final result. We present several examples to illustrate the above results and the corresponding computational savings.
 Thank you very much for coming with me today. And I'll give two lectures. The first one will be related mainly to multilayer Monte Carlo and the second one we will touch upon multi in the expressions of the things that I'm going to say today, but not in an adaptive fashion as we will see. So first go then to discuss a little bit about the activity within multi-level Monte Carlo. Let's see how this goes. So the plan is as follows. I will, well, these are some pictures from the place I come. This is a chaos actually. It's a King Abdullah University of Science and Technology. Some of the students are here today. And let's see. This is a plan. So most of you are very well acquainted with Monte Carlo and possibly multilayer Monte Carlo as well. So I will recap a little bit just to set up the notation at the beginning. So we will speak the same language. And then we go into different versions of our activity. As you would see, there are many ideas for adaptive approximations in multilayer Monte Carlo that we can see. And for that reason also, I will touch upon different models. I mean, one of them will be computations with random PDEs, that is PDEs that have random coefficients. And the second one will be diffusions. These are etostochastic differential equations. And the third one will be pure jam processes, say stochastic reaction networks, whose trajectories are piecewise constant. So I think this will motivate a little bit different ideas that we will have in this. And finally, we will discuss a little bit of the contribution multilevel and optimal hierarchies. Hopefully, if the time allows, I will also discuss optimal stopping with non-asymptotic rules within Monte Carlo. So I'll try to really be adaptive in this process. So let's begin with our first example. This is a linear elliptic PDE that is very popular in the context of random PDEs as a model case. So here, the coefficients a and f are random processes, stochastic random fields. And mainly, we'll concentrate on the effect of a, because a is with the two, right, the one that acts non-linearly on the solution u. And out of the solution u, we are interested in computing a statistic, which is real value. Again, as an example, in this lecture, I'm going to go through the computation of the expected value of a functional of the solution here depicted by psi of u. Psi could be linear or non-linear. We have some restrictions on smoothness on psi. But for the moment, think of a linear functional that is the magnetic acting on u. Again, this is not the only thing that you can compute with Monte Carlo, but this isn't the one that we will be focusing mostly through the lectures. All right. So the goal is set. We want to compute this expectation. Now this expectation involves the solution of a differential equation. Seldom, do we have the actual exact solution of u of our disposal? So instead of u, we're going to approximate it by uh, say, where uh comes through a discretization of this equation. In this case, of course, it will be only space discretization. But in general, you may have time as well or even other parameters. So think that you have prescribed the use of your favorite method. In this context, you could define it elements, find a difference, find a volume, so whatever. And that discretization is indexed by the parameter h. All right. So h will govern in two aspects of the solution. One will be the approximation characteristic of uh to u, but at the same time, the cost of evaluating uh. With the assumption, of course, of consistency, meaning that when h goes to zero, uh converges in the proper space to u. And at the same time, the cost of evaluating uh goes to infinity. Right? Otherwise, we will let h go to zero and get the exact solution. Now, since we were interested in this function actually, which is a real value random variable, we will approximate, I mean, the notation denote the approximation by gh, where now h is the discretization parameter. All right? So that's hopefully a natural notation. Let's go through the error then. The difference between what we would like to compute, which is the expected value of g minus what we actually compute, which is a sample average based on IID samples out of the inexact quantity, which is gh and not g. Right? So we split upon two main sources. Right? One, again, the trick is to somehow subtract the expected value of this approximate quantity to split the difference between this and this into two things. Right? One that we call the bias error, the other is a statistical error. And for the bias error, we have essentially an estimate of this kind, which is coming directly essentially from the deterministic analysis, right, point-wise, and then integrated upon this expectation. So essentially what we have is that if your method has order w here with respect to h, which is the discretization parameter, when essentially the constant that comes from the deterministic analysis will be random, of course, because it will be depending on a and an f. But once you take the expected value, you get the average of such a constant, which is no longer stochastic, it's just a deterministic constant. All right? So this analysis holds as long as there are no changes of order with respect to the randomness in the problem, which is a reasonable assumption. But if that happens, actually, this will not be the correct way to draw things. But again, this is an actually way to explain how this works. So let's again consider that this w doesn't change with respect to omega, and therefore we just get a constant here times h to the w. Okay, so this is more or less a standard discretization error. Second term is the standard statistical error, right? And this error essentially is mainly governed by the number of samples in. And I say mainly because the random variable that you're actually sampling depends on h, right? So this is an unbiased error, correct, because we take expected value here and then some sample of the same random variable in an ID fashion. So we only need to control its variance, right? So and that is motivated by the central limit here, essentially properly scaled, right? This will be approximated by a normal distributed random variable. And this allows to write then later the constraint of having this error smaller than a certain quantity in probability, you know, as a constraint in the variance, right? So as we can see, again, we have two contributions here, one coming from h, the other coming from m. And as long as the problem is a rational well-posed problem and there's nothing strange, these variants is uniformly bounded with respect to h, right? And therefore they claim that m is essentially the quantity that governs the statistical error, okay? And if the problem has a singularity of some kind, maybe this is not completely true, but I mean, again, in the examples that we're going to visit, this variance is going to be bounded uniformly with respect to h, okay? So two sources, one is governed entirely by h and converges algebraically with respect to h. The w here is determined by the numerical method that you use and the regularity of the problem that you're solving. And this m is very robustly coming into this analysis, right? We only need to ensure that this variance is bounded, essentially. So now that we have the two main errors identified, we try to see what is the corresponding complexity, which is the computational work that is needed to achieve a certain accuracy with Monte Carlo in this context for a prescribed confidence level, okay? And to do that, essentially what you do is to minimize the work subject to the total error constraint, okay? So the total error constraint, again, improbability, because this was essentially an analyzing probability, is the sum of the bias plus the statistical error. The bias was governed by this, the statistical error governed by that. And the total work, naturally, is what? Is the product between the number of samples that you have to solve for, right? Sometimes the expected cost per sample, okay? And again, this is an example, right? And the example that we have illustrated is a linear elliptic differential equation, right? And in that case, what we have to do is to solve a linear system every time we need to create a new sample of the random variable. So if v is the dimensionality of the space that you're solving, right? I mean, could be 3D, for example, 2D, 1D. Then h to the minus d, if h is a discretization parameter, is essentially the number of degrees of freedom that you have in your discretization. So h to the minus d is a degrees of freedom. And then to the power gamma is essentially the cost that one incurs by solving the system, okay? So gamma equals one will be an optimal solver, and gamma equals three will be the Gaussian elimination, okay? So gamma can raise between one and three. And of course, h to the minus d is the number of degrees of freedom in your system. Is that clear? Okay? So this adapts to, again, to the type of solver. If you use one type of solver or another, this gamma may change, okay? Let's see what else is to be said. Yes, okay. So I mentioned the following. I said that this cost was understood in the mean sense. And why do I say that? Well, if actually you have an iterative solver, right? And the spectrum of the operator that you're solving with is random, right? Because the a is the feasibility coefficient, was random itself, the number of iterations may depend on the realization, right? So what you're really interested in is the mean cost, not just the cost. I mean, the mean cost will be representative. And there will be variations with respect to each of the realizations. Is that clear? Okay? Now, if you use a direct solver, this will not be the case. But okay, I just want to emphasize a little bit on that. Now, if you actually do the minimization of this bound and subject to this constraint, you get this classical cost with respect to the tolerance that you're imposing. Okay? In the context of the central limit theorem, the confidence level is only acting on the C2. So it really plays a role in this constant. I mean, in this exponent, it's just acting on the constant. And that's why you don't see it. But again, everything now depends on the accuracy you want to achieve and can be split into two contributions. One that is told to the minus 2 and the other is told to the minus d gamma divided by w. Now the first part, told to the minus 2, is the classical cost in terms of number of samples of a Monte Carlo method. This is essentially what is needed to make this less and tall. So that will be the cost that one will have to pay if each sample cost me over 1. The second part of the cost is coming from the fact that I need to make h sufficiently small so the bias constraint is satisfied, meaning that this term is actually less and tall. Okay? So essentially what you see here is the number of samples times the cost per sample, which is essentially what you were seeing here as well. Okay? Now, what Monte Carlo tries to do is to go beyond this bound. So what we will improve somehow is this part and this will still be a barrier that we will not be able to cross at least with Monte Carlo methods. So just to give the notation again for multilevel Monte Carlo, let's recall this construction of Henry and Giles. Let's take a constant beta for subdivision. For all purposes, think of beta equals 2 for the moment. Okay? So what we say is the following. Let's just take in a single h that was meant to satisfy the bias constraint. What we're going to do now is to take a number of such h's. In principle, if you come from analysis, from numerical analysis, this is a natural way to proceed because every time you are actually trying to solve a differential equation, you try with a certain mesh. If you're not happy, you subdivide it right and you test the stopping criterion. And if you're not happy, you subdivide again until you meet the stopping criterion. Normally, what you do is that once you get to that level that you accept the error to be smaller than your prescribed accuracy, you throw all the other things, right? You don't care, essentially, and you just work with the last accepted guy. Here on the contrary, you're going to keep those for a certain purpose. And the certain purpose is to produce an effective control variant. So what one is doing here, again, let's go through the definitions first, is to introduce these delta operators that are differences between consecutive approximations. So essentially, this, remember, is a functional evaluated using a mesh of size hl. And this is a functional using a mesh of size hl minus 1. This is the difference between two random variables. And we are assuming, of course, that they are not independent. They are going to be sampled according to the same omega. That's fundamental in this construction, because that is what will make this difference actually go to zero, essentially, when l, bit l here goes to infinity. That again comes from consistency. In l2, this is what we're going to mainly need for this to be working fine. So let's go again to the construction. The construction of these hls is done via a geometric discretization, meaning that you start with a given mesh. If you're not happy, so divided by 2. If you're not happy, so divided by 2 again, by 2 and so forth. That creates this hl based on h0 divided by beta to the minus l with beta equals 2. Is that clear? So now that these deltas have been defined, then one rewrites. Remember this was the guy that we were trying to use as an approximation to this one, so that the bias between these two was less than tall. Now this is somehow the random variable that one is going to apply control variance to, and before doing so, one rewrites this expectation as this other sum. So this is an exact representation, and the sum is running overall levels of approximation from zero to the level l capital that satisfies the bias constraint. Yes? Now instead of using multicadload here, one instead uses multicadload independently for each of the terms coming into this sum. All right? So, and that's sort of a scary proposition because the first time you see here, this is saying at the beginning I only had one expectation to compute. Now I have l capital plus one. How is it going to work? Well, okay, before we look into that, let's see how this looks. And as you see, this, which is actually the estimator for multicadload, looks precisely like the sum of l capital plus one multicadload estimators, where each of them now is sampling something different, each of them is sampling these deltas. And of course, remember that for l equals zero, actually, this is not the difference, it's just g evaluated on the course's approximation. Okay, so you can think of this as, you know, g is h0, right, if you think of hierarchy approximation, and then the other things are just corrections, right, that are being added to the course's level. Now that we have this, how do we motivate, at least first intuitively, why is this going to work better than the original multicadload approach? Well, let's see, the bias is not going to be changed. We said, essentially, if this L is chosen so that, you know, this bias is less than told, what we are rewriting, essentially, is this quantity. So if you're going to actually improve this sampling is because we're going to pay less in the samples, essentially. The difference will not come from the bias. Yes, please. So how do we choose for each level? Oh, this is coming. This is coming. This is coming, don't worry. This is coming. Okay. Now, the variance here is the variance of multicadload, right? So what did we say that the constraint, you know, the statistical constraint looked like? Well, essentially, thanks to the construction, based on independent samples, the variance of this estimator is just one over ML capital of the variance of this GL capital. Now, you think that the GL capital is more or less close to G, so more or less you can expect that this looks like the variance of G divided by ML, right? So this is essentially saying that in the deepest level of approximation, the capital, the one that satisfies the bias constraint, you need to use total to minus two samples, right? And that is what you expect from multicadload. If you, on the contrary, use multi-level multicadload, right? The variance for the estimator, again, by the construction of independence, not only between the samples that we use on each level, but among levels, right? We get some of variances across levels. The first variance, the one that is coming from the level zero, actually, looks very much like this term, right? This is not a difference. This is not a difference as well. But the other terms are actually one over M, right, times variance of differences. And what's the point here? Remember, as little L gets larger and larger, we hope to have convergence, okay, just by consistency. Now, in which sense we expect to have convergence in this setting, what is most useful is to consider not only you have convergence in the mean, but also in L2, which means that these variances of differences are going to zero, right? This is where multi-level multicadload becomes effective. And if they actually become very, very small, it means that for really deep levels when L is large, right, each sample cost me a lot, but the variance is correspondingly small. So the number M that we have to use is really small. And this is where you get real advantages, because you are sampling much less in the most costly levels, all right? Now, you have to check some conditions usually to make sure that the choice of the level zero say sound, but again, the whole thing can be put in terms of these inequalities that again indicate that L2 is very relevant in the choice of the multi-level parameters. So far so good? Any questions? You said there was to be correlation between the different levels. No, there's not. Actually, no, they are constructed totally independent. Okay? That's important. Yes. From the independence in the construction, you get immediately this sound. Okay? Otherwise, you will have extra terms, all right? It's not that it's impossible to do that, but the construction here is done on independence in some chance. Okay. So cartoon type of approximation essentially, we are somehow moving samples from sort of deep refined levels that are costly into levels that are cheaper and less refined. Okay? Okay. So if you want to go beyond that intuition, you need some assumptions, right? And these assumptions may or may not hold in your problem. You have to verify that they are satisfied and sometimes you may have precisely this kind of convergence and some other type, I mean, depending on what problems you have. So here we are again going to concentrate on this algebraic type of convergence with respect to H. Remember, H was the parameter of discretization. You can think of delta x, delta z, delta y in your linear elliptic PDE. And we were assuming that the HL was essentially H0 times beta to the minus L, right? A geometric subdivision. So this first equation here is nothing different from what I wrote before, right? It means that the error in the mean is H to the w, essentially. You see? So nothing is different. It's just a rewriting of what we have before. Now again, this w, which is the exponent of the weak convergence, depends on the assumptions on your problem and the numerical discretization that you're using and has to be found for each problem, right? The variance assumption again is indicating some kind of error rate, but this error now is measured in L2. And instead of taking a difference between g and gl, now we are taking, remember, the difference between gl and gl minus 1 when we put it inside of the difference. And we have a strong rate of convergence that is s, which is not going to be w, right? In general. So for the last one, we have the cost, the mean cost of evaluating a difference on the level L, right? So this is delta gl. And this again entails the cost of evaluating gl and gl minus 1 because we're going to take the difference. Usually when assumptions that gl is the one that dominates, so this goes like that. And this is nothing from what we said before again, because before we said that the cost was H to the minus d gamma. And again, this is the same thing. Just rewritten in terms of this beta and this L. I emphasize again that the constant c, w, and s also depend on the problem and have to be found. And remember, gamma may depend also on the solver that you use, not just on the problem and the discretization. So if you change the solver, gamma will change even though you keep the rest of the things constant. As an example, for a smooth linear PDE of the type, linear type that we discussed today, if you use multi-linear piecewise continuous finite elements, we will have 2w equals s equals 4 and of course, depending on the solver, use gamma ranging from 1 to 3. So what's the total work? And here we start getting close to your question about the m's. The total work will be found in the same way as we did with Monte Carlo. We have L capital plus 1 levels. On each level, we have to use ML samples and the cost per sample. It's the work of delta GL in a expected sense. So what do we do? Well, remember, we have to find ML. So if you fix L, essentially, you can write an optimization problem and I'm going to go into that a little bit more when I go into multi-index Monte Carlo, but we can just say it again. I mean, remember, if you fix L, the bias constraint is not going to be important and what you want to minimize is minimize with respect to the ML's from L equals 0 to L capital. What? The total work. And this is the sum from L equals 0. I hope it is possible to see this. There's no reflection or anything like that. So L capital ML and this work delta GL. And then, of course, there's a constraint. The constraint is the statistical constraint. The statistical constraint, again, is the variance being less than the total square. So what is the variance again is just 1 over M0 say variance delta, actually G0. Just write it differently. Plus the sum from L equals 1 to capital of the variance delta GL. And this is less than total square. Of course, there's a constant here in between, but I'm not writing it. The constant comes from the confidence level that you want to impose. Okay. So the control variables here are ML, capital plus 1, right? And here you see how they act. So you can just write the last of relaxation and you realize that you can find these ML's. And ML's actually are going to be, at least if you make a real value relaxation first, this will be total to the minus 2. And then there's a constant here that is coming out of the levels. And this essentially is the variance of delta GL, I think. You may correct me if I'm doing this wrong. And then this is the cost, essentially. And this is for root here, if I'm not mistaken. Okay. So you can first write this lagrange of relaxation, solve it, and then say round it above, I mean, take the ceiling of this number to make sure that everything goes well. Okay. And from there, you get to that formula, essentially, assuming that the cost is not trivial, essentially, that you're not getting less than one realization per level. But that is a technical assumption. You can stay with this for the moment. So what do you see from here? You see that after you substitute the optimal number of samples, the cost only depends on L capital. Because these are problem, say, and discretization-based constants. So once you give the discretization method and the problem, of course, this work L and the variance L, essentially, are given. Maybe you don't know them, but precisely a priority, but they're given. And this is a sum of positive terms. So first thing that you observe is that when the tolerance goes to zero, L has to go to infinity. Right? Is it clear why? You have to satisfy the bias constraint, right? So the bias constraint, essentially, was h to the w less than toll. And L is coupled directly with this h to the w. So L, essentially, goes to infinity like a logarithm of toll, absolute value, say, that's a ballpark. You may have log log terms on top of that, but really it's log toll. So this is going to infinity as toll goes to zero. So this sum can only do two things, either stay bounded as toll goes to zero or diverge. If it is diverging, it's because this work is somehow growing faster than what the variance is being able to converge, right? But if it is actually finite, which is essentially the case where this variance is dominating the work somehow, we have a very interesting case. But it's sort of not intuitive the first time you see it as well. It's the fact that you get a cost that is essentially all of toll to the minus two. And if you remember Monte Carlo, when we discussed it first, toll to the minus two was the optimal cost with Monte Carlo when each of the samples, essentially, in the mean cost me all one. So essentially, this means that you don't feel the cost of discretization even though you're refining the levels, right, infinitely much as toll goes to zero. Is that okay? In principle, not, right? So somehow it means that even though we still have to sample of those at the level L capital, we only sample so few, so little samples so that that cost really does not influence the rest of the cost that come from the other levels. Because we still have to sample in that way. That is clear, right? Okay, so now that you see what is coming, let's see the result with these assumptions in mind. L capital, as I told you, is determined by the bias constraint and essentially is linked to toll by a log factor, right, where the log essentially comes connected with a weak rate of convergence and the logarithm of the subdivision multiplicative constant beta. And the samples ML are chosen optimally like we discussed before just by solving that Lagrangian relaxation problem. And after substitution, essentially we get this fantastic case, toll to the minus two, where S is larger than the gamma, right? So essentially both the work and the variance are behaving with respect to L in an exponential fashion, right? So the work is growing with the D gamma rate and the variance is decreasing with the S rate. So as long as S is larger than the gamma, this sum is convergent. And this is what explains the first result. The second is a limit case where essentially they outmaneuver each other. So essentially they balance each other. You get something that is constant here. So the total result here is L on this sum, L times a constant, and then you get L squared and L was log toll. So that explains the second line here. And in the last one, well, the work is growing faster than the variance is converging and instead of getting something that is convergent, you get some algebraic deterioration. But still this algebraic deterioration is better than what you would get with Monte Carlo. So you're decreasing from D gamma here in the explosion to D gamma minus S. OK? So you're still exploiting a little bit of the strong convergence. I mean, the convergence in it too that one assumed. So OK? So far for the introduction with the uniform methods. Now the question that we want to more or less touch upon is how to extend this Monte Carlo into non-uniform adaptive discretization settings and whether it's worth it. OK? And there are several ideas to introduce non-uniform discretizations or stochastic levels, if you wish. So let's begin with the... OK, well, this is... I put just to include in this slides and I already gave up trying because, I mean, the pages are not enough actually to describe so many things. And I rather just direct people to look into mics. Himongu said to keep up with the growth of the community. And this is actually what you get the best references in the... List of references in the area, so I'm not going to talk too much into that. So let's go to the activity then. And in this context, again, we're trying to approximate the expected value of an output quantity. In this case, I'm switching the example. I'm going to stochastic differential equations, itto differential equations. So in this setting, each of the paths evolved in finite dimensions with respect to time, right? So you can do this with the sum of the paths to a final value and then in that final value, you may evaluate the function g here that comes in this example, like a put option in finance and then you have to integrate this g versus the PDF of the final time. OK, so you can do this with Monte Carlo just by sampling many, many paths and again averaging as we discussed before. So here are some references again. I'm not going to stop and I'm going directly into the problem formulation. Here is an SE. This SE, if you have not seen it before, it's not a problem because you can think of this discretization via Euler-Marujama that is coming just in a second. So x is something that is finite dimensional now, it's not a function as in the first example. It evolves in RD. And g is a given function so that our goal is to compute the expected value of g of xt and so the evaluation at the final time of the solution of an SE. Again we are given an accuracy constraint tall and a prescribed confidence level and everything is driven by k-dimensional venerable process. All right, this has plenty of applications that I will not stop upon and what we're going to do is through this adaptive Monte Carlo and actually see that there are several versions now of adaptivity. The first one we did was with respect to the weak error but later we did something that controls the strong error and yeah, let's go into that later. So first discretization algorithm. This is just Euler-Marujama so if you have not seen SEs before, this is a new term that comes in. So remember it's coming through these. If you just think of the first part, it's just forward Euler right based on this function, the drift function a and the second is just the same type of discretization. You see you have to evaluate b on the left point, essentially the current point x bar but you multiply times these random variables. These random variables are normal distributed with mean zero and the variance is just the delta t that you're taking. So the delta t does not need to be uniform at all but yeah and again weekly you can just scale them to normal zero one random variables and the multicolline estimator looks exactly as we discussed before. It's an average of IID samples but the IID samples are not from the exact process but from the solution of coming out from this forward Euler. Forward Euler of course is not the only discretization that one can do in this context but it's the one that I'm going to exemplify. Okay so as we discussed before this can be split into tolerance for time discretization and tolerance for statistical error and this is a particular case as we discussed before the cost per sample here is essentially proportional to the number of time steps. The number of time steps essentially is proportional to one over a ton. This coming from the weak convergence rate one of Euler which is the same as Euler-Maluchin in this setting and this brings us to the complexity told to the minus three for Monte Carlo. Okay okay and now you want to introduce adaptivity. So what is adaptivity in this context and this context adaptivity means the ability to choose the time steps non-uniformly and this non-uniformly can depend also on omega. So the time steps such as t0 is 0 and the final time step is t capital are not going to be equally spaces and they may adapt to omega because the dynamics of the solution may depend on omega in a dramatic way maybe. So in this context then we will create realizations that depend on x bar, I mean here denoted by x bar and depend on this mentions delta t of omega. So again why do we do this? Well sometimes there are non-smoothness on A or B that ask to, I mean to deteriorate the usual behavior of the uniform method so ask essentially for the adaptivity in a natural way and what do we need to be able to produce such delta t? Well you need, at least with our approach, you need to rewrite the weak error in this other way essentially, you need to introduce an error density which is a function rho w that does essentially not depend on delta t so that the error, the weak error can be split essentially in contributions that are coming out of each of the different discretization steps. You see here if you discretize over the time steps, this is a sum over time steps of essentially delta t n square times rho w. Now the fact that delta t here is inside of the expectation is allowing for the choice of stochastic time steps because rho w may still vary dramatically with respect to omega, that's the point. If this rho w will not be varying dramatically with respect to omega then it's cheaper actually to take a non-uniform but deterministic discretization with respect to delta t, is that clear? There's another head to do this fully stochastic so whenever possible you avoid that. Now the funny thing is that in this context there's not only a weak density but there's also a strong density. You can actually make an approximation and you see here the power is one but here the power is two, you can actually write the L2 error here in the same fashion for Euler where now this rho s is not exactly this rho w, this is still something that you can approximate with computable quantities and it's able to distribute the strong error across the different intervals that come into discretization. In principle this is not obvious because when you square this you may get a lot of cross terms. So to be able to tell that actually the diagonal terms are the one dominant this needs some understanding. Okay, but what is the whole point now is that once you actually produce these error densities you can produce single-level adaptive algorithms. So single-level adaptive algorithms will try to find these time steps so that say for example a bias error is satisfied or even say this type of L2 error is satisfied. And what does it allow then? It allows for the following construction. If you want to create a multi-level adaptive method what you derive essentially the discretization with is with a sequence of tolerances. So if you have a reliable single level adaptive method the only step that you need to go from there to a multi-level adaptive method is just to create a sequence of tolerances and to say that the first level will be driven by what will be defined by discretization of my single level adaptive method so that they reach that accuracy. The next level is again refinements on those that achieve the next accuracy constraint and so forth for all the levels. So you see now each of the levels does not have a prescribed H. This notion is gone. Each of the realizations may use a completely different mesh that adapts to the characteristics of the omega that you have in mind. But still one will be able to extract the good properties of the variance convergence, the mean convergence and the rates of the explosion of the work with respect to the level. Okay? So what is optimal at least to show results on the algorithm in this context is to refine according and to stop according to these simple rules. Essentially what you do is to compute solutions and error indicators are n for each time step. And the error indicators are n in the context of this error density defined essentially as rho n times delta tn squared. I mean and rho n could be the weak error if you're controlling the weak error of the strong, if you're using the strong approximation that's not important. But again now this is defining the contribution from times the n in a given realization to the total work. And what we do to stop the discretization is to verify whether the maximum of these guys is or not above this. And this you see is the tolerance that we prescribe for the bias constraint which is what we want to achieve divided by the expected value of a number of time steps. So this essentially is the mean error that we would like to achieve on each interval. That turns to be optimal essentially to equidistribute the error along intervals. One caveat though, this n is not known a priori, right? Because this is the number of steps that you need to solve. So these have to be found on the fly. But this is a constant and this is another constant that are essentially probably independent. So what you have to remember here is that this is going to be not the same as this C of refinement. The refinement criterion is a little bit more stringent than the acceptance. So essentially if you see that you're not going to stop because some of them are going to be larger than that, well then you are going to refine many more than just those that not satisfy this constraint. Now another technical point in this construction is the following, if you have a refinement and has not been accepted and you mark some intervals in time for refinement, you need to produce essentially the corresponding interpolation values for the delta W's that drive your discretization. And these, you know, not to introduce any biases, have to be consistent with the points that you already sampled. And this is done according to Branian Bridges, right? That is the correct way to do this discretization. So this plot essentially shows that, I mean, that for a given level, for example, you use the blue discretization, right? And this entails evaluating W at all those points. But now you need to evaluate W in points in between. And of course, they are not going to be linear interpolants of the previous values, but they will differ upon. And this difference is actually given by the Branian Bridges, okay? So it's just conditional Gaussian sampling. That's nothing more. Okay, any questions? No? Okay, so what was the main idea again? Let the tolerance in the delta W algorithm define the hierarchy. That's essentially what you have to remember. And this hopes, this works in any problem. If you have a sound error density of PEs or EZEs, you can just use this approach and it will work. Like I said before, adaptivity is useful when you have issues with A or B or the process itself may have issues with respect to, say, boundaries where it is being, for example, stopped or reflected that those actually ask for adaptivity. You may have stiffness as well, and this is also nicely treated by adaptive algorithms. So again, if you have a stiff drift, you may choose at least delta T. If you don't want to go through this posterior global error control, if you just want to drive things by stability, you just may take delta T to make sure that you remain within the stability region of the method. And that again, if the stiffness is just coming from A, it can be done quite easily. So what is interesting here? The point is if only when you do a posterior error estimate, you show things about the discretization being efficient or not in the sense that the estimate for the error essentially gets close to the error you want to impose. But in this context of error densities, you can show actually, you show much more, you can show actually that the algorithm stops, that asymptotically you have normality, asymptotically you have accuracy. And again, what is asymptotic in this case is when tolerance goes to zero. There is no notion of H equal to zero because H is gone. So the only parameter that drives the discretization here is the accuracy constraint, and the accuracy is given by toll. So when you look at asymptotics in this context, it's when toll goes to zero. And of course, we have kept the confidence level constant in this discussion. We are not letting the confidence goes to one. So asymptotic normality, you also have asymptotic accuracy. I will just tell you later what this means. But it really means that when tolerance goes to zero, essentially the error is bounded by toll. And the complexity essentially is up to essentially a constant, you get the optimal complexity. So for the ethos to classic differential equations, this one discussion that I have not made, and I just put the warning here, it's actually working fine in the theory, but it's not obvious. The discretizations that you get are not adapted, I mean, they are adaptive, but they are not adapted. Adapted means that they are not adapted to the natural filtration of the problem, which is given by the linear process. Why? It's because the density, and I have not told you maybe, no, these densities essentially involve variations of the final value with respect to the current value. So essentially, they may need, for example, derivatives of respect to the position Xn of g of xt capital bar. This quantity is not adapted. It's looking to the future from, say, the current point n to final time. So in principle, you can discuss whether these things are going to converge to the correct solution. They do. But I mean, this is not true for any adaptive method that you can come up with. And no one has to be worried, I mean, at least careful about it. Okay, so that's why. So if you do examples here, for example, you have a drift singularity, this is just going to be simple in the case, in the sense that the stochasticity is really not relevant. I mean, here we have, say, deterministic time for the singularity, and the mesh is going to refine automatically around the value of alpha. We are not telling the algorithm that alpha is there, but it does find it. You see the mesh here is refined from the level, essentially, 10 to the 0 to 10 to the minus 10, just to control the error. Okay. But again, like I told you before, this is sort of the deterministic type of singularity. So in this case, it doesn't pay too much to let delta w to depend on omega. Delta w is not uniform, but it really doesn't depend on omega. So you can split the construction here in a first phase, where you actually find all these meshes that you're going to use in all the levels. And then you just run in the second phase, material multicolor increasing the number of samples that you need on each level until you hit the variance constraint. That is effective. Okay. And actually, what you can show is that the cost of the first phase to find actually the optimal discretization in Monte Carlo is much cheaper than the actual cost to solve the problem. So that's again, actually, when you have deterministic type of singularities that you can find the mesh very quickly. Okay. So in this case, you see an advantage from uniform to adaptive. You don't see much more of a difference between uniform Monte Carlo and adaptive Monte Carlo, because essentially what happens is that this singularity is not strong enough to modify the strong rate. It does modify the weak rate. But if you remember what I said before, the complexity in Monte Carlo is mainly governed by the strong rate. I mean, there's this independence of the weak, of course, but in the inequalities that I showed you before with these conversions of the sum, the strong rate was the one that was opposing the weak, the work growth, not the weak. Okay. So, but if you actually go to stronger singularities, so now we let P actually be beyond one half and, of course, you have to make sense that this is actually has a unique solution, blah, blah, blah, but it does actually, as long as P, I think is less than one strictly. So you can see here for P equals 0.5, that was the previous example, you don't see much of an advantage between multi-level Monte Carlo uniform and adaptive, but as you start raising the P, you clearly see the advantage between adaptivity and non-adaptivity. And again, I have to tell you one more thing, this has been done using the strong density, not the weak density. So this is the latest work that we have for adaptivity. It's actually on the strong, I mean, you use the strong density to rewrite this L2 error and you now create the levels according to this strong density. That's how this plot is being created. One more thing, okay, the alpha here was uniformly distributed between one fourth and three fourths, so actually the measures have to be stochastic, because the position of the singularities is being just jumping around. Okay, so I told you essentially how to do adaptive multi-level Monte Carlo approximation, yes please. Sorry, alpha is independent of w, that's right. So essentially you set alpha, you run the things, alpha is independent of w, correct. Yes, I mean, this was just a simple example. So people can actually reproduce and see what is going on. But yes, I mean, things can be more complicated because the singularity, as you're saying, can be coupled with the dynamics, okay. So it could be that actually the singularities appear because of some drift problem, right, with the stiffness, or it could be that x is getting close to a boundary and the boundary has an issue. So you can think of the boundary being penalized and penalization is drift stiffness, right, so yes. But that will not be able to be written in one line essentially, so yes, at least this motivates what you can gain. Thank you. Okay, any more questions at this point? I'd like to then summarize what we did on this topic. So I hope the idea of adaptivity in this context that I gave you is quite general. You can use it for any problems you want essentially, as long as you have a sound single level analysis, if I may say, okay. The message is that it is better in the context of multilevel to use an analysis that will go towards the L2 error, if possible, right, drive the analysis of the construction of the levels according to that. Of course, you need the bias error to control the number of levels you're going to use, but at least the distribution of the shape of the meshes should be done according to the strong error. And actually that can be proved, actually, that is optimal to do so. You just write a problem that is related to these, but instead of writing this part here, you just write the error density times the, I mean, the integers of the error density times the delta t's, you just write the whole Lagrangian relaxation, which now will involve finding functions, functions delta tL, say, from L equals 0 to L capital, and then you will see that delta tL is directly expressed in terms of the row strong, not the row weak. The only, the last level is connected to the row weak, but all the rest are driven by the row strong, and therefore that motivates directly to the use of the strong for the activity that we found later, sorry, but, I mean, the first works we did was with the weak activity. So how do we do with the time? We split in between or? As you like. Okay, I can just draw the conclusions and it would take five minutes. Let me put the conclusions and then, okay. So as conclusions, I told you again that this way of thinking for multi-hebel and adaptivity is quite general, you can use it anywhere you want. You have to be careful, though, that when you replace the statistical constraint, which means that the statistical error has to be less than tall with a prescribed probability, computationally it's much more effective to write that into a variance constraint, and this can be done as long as you have a central limit type of approximation. The central limit does not apply directly in the multi-hebel-vontical setting because you have some of independent but not identically distributed random variables, so you have to use resort to some linear or failure type of result, and you have to verify again these four each of the problems that you work with. So you have to be careful with that. Okay, so let's leave the puzzle random measures for the next one. Thank you. Five minutes? Yes, sir. Five minutes. Thank you. So we are going to discuss another idea for adaptivity in this second context. The second context involves paths that are also random and depend on time, but as opposed to the stochastic differential equations that we were given before that were etotype, these paths actually can be solved exactly. So these paths are given by piecewise constant functions, and the only random thing about them is essentially the time where they jump, and which kind of jumps we have them. These processes are very popular, they are macro processes of course, and they are used in Xenium applications, but let's go a little bit more into them. This is a logarithmic scale, this is time, and here in blue, red, and green are three components of this vector that is changing according to time. The numbers that are shown here involving are the number of species for each of these three components, so you can think of three populations randomly changing inside of a system according to maybe biological reactions or chemistry reactions for instance. Even though these green paths do not show it, they are as discontinuous as the blue ones essentially. It's just that you have logarithmic scale and you may be fooled by the logarithmic scale, that's it. So that's a motivational example, this is a gene transcription and translation you have here in this system, the component G that may react given G itself an M, M that actually may react again and give M and P, and the P that actually goes into D when it meets itself. So you see that there are sort of reaction channels, essentially when this thing happens, you know deterministically what is going to happen. The only problem is that the time where this occurs is not going to be prescribed, a theory is going to be random. And let's think of a simple example, let's think of water say in a container and water can be created from hydrogen and oxygen and it could also dissociate into hydrogen and oxygen, so this is a reversible equation. But the way the speed of these reactions actually happen depends of course on the temperature, but if you fix the temperature out of this discussion, they depend for example on how much hydrogen and oxygen is there because if there is no hydrogen, oxygen will never meet hydrogen to create water, right? So the likelihood of these events to happen depends on how dense these components are in the system. So and this can be mathematically described. So what we want to do again will be to compute some given output quantity of interest at a given final time up to a certain tolerance again with high probability. So how would describe the evolution of such a system? I told you before this was an X that lives on a lattice essentially because it's counting. So the components of X essentially describe different species in the system. So we count how much of the species one we have, how much of the species two we have and so forth, okay? So X is taking values on these lattice, right? And it was according to time and of course omega. Now we assume that we only have finitely many possible reaction channels like in the case of water for example, we may have three components, one is how much hydrogen, how much oxygen, how much water, right? And we can only have water dissociating into hydrogen and oxygen and the other two meeting to get water. So when this oxygen and hydrogen meet, right, we assume that they can only create water. So there's nothing stochastic about the fact when we know that they meet. The question is when they meet and this is the stochastic part. Okay, so when this thing happens, right, the J reaction type of in the system happens, the only thing that we have is a change from X that was the state of the system into X plus nu J. In the case of this water for example, when oxygen and hydrogen meet to create water, you have a plus one in the water and the corresponding minuses in the oxygen and the hydrogen. Now, when does this happen? Is given by the so-called intensity functions for the jumps or in this context also called propensity functions. These are non-negative functions, A, J, that are evaluated on the state X and give values on R plus. So that the probability of having a reaction of type J in an infinitesimal interval that goes from T to T plus DT, given that you have the state X is essentially A, J of X times DT. So for instance, if you have no hydrogen in the system, you know that new water is not going to be created. So this A, J corresponding to that, that reaction will be zero in that interval until that situation changes. So again, we want to compute the expected value of a given quantity for the given of several G and the applications are everywhere. So I'm not going to stop on that. So we have plenty of references and as well this is not complete, but we just wanted to show a little bit the efforts on this multilevel version. So what is interesting about these systems? First thing as opposed to the previous discussion that we have is that exact algorithms are available. What do you mean by exact algorithms in this context? You remember when we sample a la Monte Carlo to approximate quantities of expected type, we always have statistical error unless the problem is trivial and has no variance and then we are not discussing anything stochastic. So exact in this context means pathways exact, essentially that with finite work you can simulate a process that has the correct law. So there is no bias error. You can do unbiased simulation with all one cost, which means that already you have the beautiful regime that we discussed before that was total to minus 2. The cost is total to minus 2, the number of samples in Monte Carlo and then all one to solve for each path in the mean. So the whole discussion here is about constants. So you do variance reduction in order to reduce the constant in front of the total to minus 2. So that's the whole point in this context. As opposed to the other case where we actually can push the complexity from total to minus 2 log square times total to minus 2 into total to minus 2. Here it goes from total to minus 2 into total to minus 2. And still you can get implicit gain. So constants matter, that's what I'm trying to tell you. So then if you have an exact algorithm, why are you going to introduce an inexact one? Well, because the two things, because the statistical error is already there. So having a bias that is zero in the presence of a statistical error may not be optimal because the cost of simulating such exact path may be a lot. So introducing a path that is a little bit cheaper and has some bias and the bias is smaller than the statistical error will not change things dramatically from the error point of view, but it will do it from the work. So it is actually a good idea to actually introduce bias in a bias simulation. Right? Good. So what is the problem then here is that essentially the cost comes from all these events that are happening in the system. So if you actually have a lot of molecules in a big reactor, you can think that there are a lot of reactions happening. And if you have to track all of them to do this exact simulation, the cost is going to be horrible. Okay? So, well, essentially these interarrival times between jumps are exponentially distributed according to, with intensity, there is a sum of all the AJs in the system. And this can be a large number, depending on the system. Okay? So as opposed to that, then one would like to do approximate time stepping. Time stepping that is cheaper but introduces a bias. This has been done in, I mean, pretty much from the tau-lib, there is a forward-oiler method to more refined versions of it. And you have drawbacks, of course. First is you're going to introduce a time discretization error that you didn't have. Right? And two, you may have non-intuitive, at least qualitatively wrong type of behavior. Essentially, if you go back to the discussion that we had at the beginning, our process is actually counting. Counting molecules of different types in a system. These are not going to be negative, right? I mean, if you're doing the proper thing, the true trajectories will never cross zero. But if you do these sort of approximate algorithms with a non-zero probability, you may cross these things. And you may want to avoid these events because they have an issue on the accuracy. The global accuracy that you can get. So there are ways to deal with these. For example, the pre-lib is just to adjust adaptively the time step, you know, the control, the one-step exit probability. But we would like to always relate things with respect to the goal, which is to relate the old type of errors that we have in this approximation with respect to the goal, which is the expected value of g of xt capital. So what we did in this line of adaptivity was to first create a single-level hybrid algorithm that at each time step, adaptively switches between the exact path simulation and this tau-lib sort of oiler in order to do what? To minimize the computational work. So why is this somehow needed? When you're close to this zero level that you don't want to cross somehow, the tau-lib, which is the approximate forward-lowered type of discretization, does have a cost that is not being reduced somehow. So you're not gaining too much from using it. But the exact simulation essentially has a cost that I told you was essentially going up with the number of particles in the system. If you're going to get closer and closer to zero, that cost actually gets smaller. So close to boundaries in this respect, the exact simulation is advantageous because it doesn't really use bias and it's cost to do because the inter-rival times are relatively large at least. So this can be quantified and a rule that switches between these two things can be done. So again, this is adaptivity, but now not adaptivity just with respect to time, but between methods. So you're combining two methods. One that has sort of infinite accuracy or infinite order with respect to time discretization and the other that has order one, which is essentially coming from Euler. So that was one thing. The other thing is that this idea can be generalized into multilevel Monte Carlo algorithms for our control again. And not only that, it's like you can do this according to the different channels. So I told you before that when you introduce this, there are different J channels. And with respect to these Js, what matters really is the propensity. If the propensity is sort of large, you're taking a lot of events over a unit time step. If the propensity is slow, these things fire not often. So if you can actually split among the reactions, do a time splitting method and essentially treat some of the reactions with the approximate method and some of the reactions with the exact method. This can be done as well, even in the context of multilevel. And again, all this, what is going on here is that you're introducing biases and approximations that have much, much less cost when compared to the unbiased algorithm. This is why this gets useful. So how do we go about it? To motivate the different methods, the whole thing can be written through this course random time representation as accounting algorithm. Essentially the value, the number of species that you will have at time t, essentially is the number of species you have at time zero, right? Changed by what? Changed by the number of times the Js reaction occurred. And the number of times the Js reaction occurred is described here by a Poisson process whose intensity is non-trivial, it's just the integral from zero to t of A J, the propensity of the reaction J. This is again a non-negative function and this is just a random time change on a given Poisson process. So this is a natural number that goes from zero to whatever and this is just the vector that introduces the changes according to that reaction. Is that clear? So it's just bookkeeping, you see? And if this A would not depend on X, this would be trivial to simulate, right? You will just change accordingly the time in a unit of Poisson process and know exactly when the jumps are going to occur in your system and that's it. So to simplify this, when you want to do approximations, you just do a donor method, essentially you freeze the value of A at the left point and approximate this increment by a Poisson random variable with intensity A J times delta t, if you wish. And the first time I wrote it, they put tau here and this is why they call it the tau-leap method. Okay, so first problem that you see here from this approximation point of view is to find tau so that you're not going to cross these important boundaries with a given probability. Delta is given. So this tau will depend on X and delta and this will prevent you from jumping according to this inequality. What we did actually is to do some large deviation, turnoff type of bound to actually come with a computable formula that one can use to ensure this inequality. And then from this delta t, we can actually discuss whether it is worth taking the approximate step or the exact one. So and then one step farther is to go into a global probability of going out because remember this discussion that I just made was based on one step, essentially. If you're here, how you have to take delta t, so the probability of exiting in this time is below a certain value. Now what you really want is the probability of the whole path going out or not going out. So you have to link the one step probability to the whole path probabilities and then link these quantities into the effect of the global error that we really want to compute and approximate. So this is the first step that is good because I mean in principle the probability of going out for a path, if things are going to be done correctly, it's a rare event. So you want to, instead of sampling that rare event, you want to compute quantities that are sort of easy to deal with. And this approach does it actually. Okay, and like I told you before, then one can write a multilevel Monte Carlo based on these type of paths, right? Find the m's and find the l. So you control the exit type of error, essentially that is related to the indicator function of the path going out. That will produce an O1 error. So essentially that's why you have to compute to control this via the exit probability. Then the weak error that is more or less what we have discussed before with the diffusions and then the statistical error that doesn't change from the discussion that I made before. So the issue of this is here is that you have to do things correctly to simulate couple pairs that have the right law and are correlated in the right way. And of course being able to accurately estimate the different quantities in the multilevel construction. This is particularly nonstandard here from respect to the diffusions that I discussed before for the following reason. The paths, right, they're taking values on a lattice. So when things are going to converge, right, I mean as you discretize the L, the delta t or increase the L if you wish, the paths are going to be closer and closer. But you know, in diffusions closer and closer may mean you're in a continuous space, so they may mean actually getting closer and closer. But here you have a lattice. So what does it mean to get close in a lattice? You're far or you're there? So actually closer in this setting, it means that the probability of going far is going to zero, right, and the probability of being equal is getting larger and larger. What's the problem with that? That the problem is that this creates immediately high kurtosis. And high kurtosis is bad news in the multilevel Monte Carlo context because it means the estimate of variance is difficult. So one thing that appears here naturally as opposed to the diffusion case that I discussed before is you have to work a little bit farther to estimate the right variances in your system. Remember the m's are driven by the variances. So you need good estimates of these quantities. Otherwise your estimates on the m's are wrong. And if the estimates of the m's are wrong, bye-bye Monte Carlo, right? So that's important from the practical point of perspective. I think it's in delta and the index and the... Little n, yes? It's in delta that you use to balance the probability of s. Yes, everything goes there. So you know that. Everything goes into there. I mean, each level in principle here contains a delta Tl but contains also a delta L. Because remember, there was a... There was this probability of going out, right? So you may prescribe this delta L as a function of the L and also a delta L maximum per level. So this thing goes as before like delta T0 divided by 2 to the L, for example. But this doesn't go that fast. This is chosen in a different way. So these two things combine to give the method. Is that clear? So these two are varying. All right? Thank you. Good. Thank you. Okay. So this is an example, the example that I showed at the beginning. And here you see that from the exact method, essentially, which is SSA, the error bound versus the tolerance essentially goes like total to minus 2. So essentially, this is the complexity that we have. And if you use this hybrid multilevel, you can still keep the same complexity and what you gain is a constant but the constant can be substantial. All right? In this case, it's 10 times faster. And this is actually the usual plot that you would see in our papers. It's the control of the weak error. So what the user gives in this context is a tolerance but he also prescribes the confidence level. So in this case, the confidence level was prescribed to be 5%. So when you do your simulations, essentially, you expect to see around 5% of them going above the prescribed tolerance. So you see here, the tolerance is getting lower and the corresponding sample error for the different runs of the material, Monte Carlo, distributed like this and then essentially these numbers here are the percentage of those above the bound. Okay? So you change the quantity here, these exceedance values should change as well. All right? So you can actually have, depending on the problem, different type of gains. I'm not going to discuss this too much but you can go above 10, of course. You can get 10 to the 4, things like that. But the plots look like this always. You're improving the constant. You're not improving the rate. The rate was optimal to begin with. Right? So it is actually a pain to introduce a bias approximation with respect to the unbiased method. Okay? And here's an example of our reaction splitting method where you actually do this differently in the different reactions to get even farther. Okay? So this is, I think that's it. So summarizing here, this is a different type of idea for multilevel adaptive, right? You're switching between methods. Okay? You can also do the delta T adaptive as well but here you're switching between methods. And you're going from an unbiased method into a bias method to improve the constant in the complexity. All right? Any questions? Okay. So now we discuss discretization from another point of view. We're going back to uniform discretization but now we ask a different question. We say we have been driving our discretizations with this type of geometric type of refinements, right? I mean you get an initial discretization parameter, then you divide it by a given number, then you divide it again by the same number, and so forth. Okay? This is natural within, again, numerical discretization, right? You try with a mesh, then you start to divide it. Try with another mesh. It's not optimal. It's not optimal. That's a natural question to ask, right? I mean, should we choose the HLs or the delta TLs like this or should we not? Answer? No. It's not optimal. Okay? So with some assumptions you can show that this is not optimal. This is the work we did in this optimization of hierarchies for Matias Monte Carlo. And this other work that is here, this continuation of Matias Monte Carlo addresses this problem, right? That these variances that you need actually to drive and get the optimal ms are not known a priori and somehow have to be found on the fly, right? So you need to create an algorithm that not only is adapting the ms but is trying to find these parameters that are needed to find ms, okay? So that's what the continuation is doing. And this is implemented in the lecture in this code, MIMC-LIVE, this library, MIMC-LIVE. So let's go to the assumptions then to focus. So what do we have here? We have a Matias Monte Carlo method as before, but now the levels, GLs, are going to be driven by these HL parameters of discretization and we do not say that they are geometrically distributed a priori, okay? So we're going to let them vary. So we're going to take L capital plus one of those and we are going to assume that essentially the weak rate is algebraically convergent. The strong rate again is algebraically convergent. And the work, for example, is algebraically divergent. That's in the linear PDE that we have before. And then what we have to find are not just the ms but also the Hs, okay? Is that clear? So as you know, the constants and the Qs here will depend on the problem and the solver, as we discussed before. And these are just a couple of examples, right? In Euler-Maruyama you have Q1, Q2 equals one. And if you have the other PDE example, we already discussed that the Q1 was two and the Q2 was four. Gamma again is defined by the solver and it's obvious in Euler-Maruyama, it's just one. So I don't stop into that. I just try to show you that these things can be found. And if you have no theory at all, right, if you're desperate, you can just try to fit something to the convergence as you see, right? And or the explosion in the solver. So you can actually measure things as well. So two useful quantities here is the quotient between the strong convergence and the explosion of the work that is chi and the eta is the quotient between the weak convergence rate and the explosion of the work. Of course, the best case for multivariate Monte Carlo is when chi is larger than one, the limit case is chi is equal to one. And the sort of algebraic improvement with respect to Monte Carlo is chi less than one. All right? And again, you see that if you change the solver, you may change from chi greater than one into chi less than one as I described here. In the PDE example, if you use an iterative solver, you get 1.34 and if you go into the direct solver, you get 0.89 already. So you can influence the convergence rate of multivariate just by changing the solver like that. Okay. So the optimization, how does it look? It looks like this work that was the work that you discussed before, subject to this constraint. This is a bias constraint less than this fraction of the tolerance. And this is the statistical constraint, which means that the total variance, right, is less than theta times tau divided by c alpha. And c alpha, again, is the confidence constant that is coming from the user, essentially, when he imposes the confidence level alpha. Okay. So the good news is that you can solve this problem. And of course, you can even verify the normality, how far you are or how close you are when you do actual computations. And here's a result somehow for chi equals one. You actually have, is a case where you actually have geometric discretizations. But I have to warn you, these geometric discretizations, you see, h a is somehow beta to the capital minus l. So it's a subdivision that we discussed before. This doesn't depend on the levels, but it depends on the total. So it's not your usual geometric discretization that you would think of. Okay. And the beta, actually, is not two or three. I mean, it's something that may even have an influence from total. Okay. So the spacing between the levels is not trivial either. All right. This is one more thing that is interesting here. The optimal choice of the splitting parameter, which is the fraction that you devote to the bias versus the fraction that you devote to the statistical error, actually goes to zero, sorry, goes to, oh, come on, goes to, I mean, the fraction of the bias goes to zero as a total goes to zero. Yes, that's right. So essentially, if you see here, theta is the fraction that you give to the statistical error and one minus theta to the bias. So theta essentially goes to one, why, because theta is fixed and L goes to infinity. So this theta is just going to one. All right. So asymptotically, you're putting most of your error, you know, your error allowance to the statistical error, not to the bias. And of course, the level behaves like log total and with this constant. But for non-Ka equals one, things become less simple. You can still find the relations on the theta and you can find bounds on the number of levels with respect to the different parameters as this shows. Okay. Nothing really changes from the classical results in terms of exponents, but you can identify how the input constants affect the actual complexity in the actual constants. So you can play tricks to reduce these constants knowing that they exist and the shape they have. Okay. So the shape they have in these constants in the different cases are like this. So you can actually see that in some of them, for example, the variance of the level zero plays a role, but in some others they don't. So here in this setting, if you're here, you may get an advantage by doing a special trick for level zero. So whereas in the other cases, you may not. Okay. So like I said, the exponents are the same as you have in the classical theory, but now you identify the constants as well. And you identify how you have to space the levels too. And that is more or less the discussion that I have here. These are not going to be necessarily nested because you have this beta that behaves in a strange way. And the H0, also the first mesh, also behaves in an interesting fashion too. Here's a discussion on the continuation of multi-level Monte Carlo. And again, the plots always look like I discussed before. You have a direct solver here and an interactive solver. You have the prescribed error. And here's it versus the tolerance. And you're checking essentially that when you run over and over and over the multi-level approximation, you get a fraction of those that go above in accordance to the more or less the confidence level that you prescribed. Now coming back to this question of finding the variances and other relevant parameters to do optimization, in the continuation of multi-level Monte Carlo, what we do is to solve a sequence of problems until you hit the problem that you want. So essentially, what happens? We know that in the best possible case, the explosion of the work in multi-level Monte Carlo is told to be minus 2, right? Best possible case. It could be even worse. And you can have a long talk here. What that mean? The worst is like this. So if you actually solve now a sequence of problems where you first begin with told zero you'll say, and then you go times told zero divided by a certain fraction, right? Call it F, and then keep going, right? Until you hit the actual tolerance you want, you can sum these works here, and you can get a choice on this fraction so that the actual sum of these sort of preliminary works does not create a huge overhead on the last one that you want to solve. That's the first idea. And then what happens here? That when you start solving these sort of axiomary problems, you start learning the constants that you need, okay? So somehow by creating these fractions, you do not overshadow the actual work that you want to do optimally, but you start sort of uncovering a veil over the variances that you need and the other things that you may need to find the optimal elements. That's essentially what is behind the continuation of Monte Carlo, sort of a motto with respect to a tall type of idea. Not really complicated, but it has to work quite nicely. Of course, on top of that, there are tricks to actually improve a little bit the variance estimates on the deeper levels, right? What happens statistically is that in the course levels you have a lot of samples, just by construction, right? And Monte Carlo produces a lot of samples in course levels and very few in the deep ones. But the deep ones are the really costly ones. So if you make a mistake in the variance there, you make a mistake on the M and things go really rotten. So to get a better hold on the deep levels, you actually do some kind of Bayesian type of estimate where you use as a prior model an extrapolation that is coming out of the course levels and you enrich it with a few samples that you have from the deep levels. That's what is behind too. So it's not just one idea, but I mean, I'm just telling you the couple that I consider most relevant. So the good news again is that if you do all this sort of little learning process, the complexity that you get is not much of an overhead with respect to the optimal line here. That will mean that for this sort of pink magenta type of a curve, you are just writing the cost that you would incur if you would know all these nice constants beforehand with respect to finding them on the fly. So it's not too much that you have to pay, really. So that's with respect to the optimal splitting theta and it's not really much relevant. So what is the second part of conclusions? We didn't discuss much of normality, but we have a result using the linear failure in this context. Better splitting between Bayesian statistical error can improve the work a little bit and this continuation of Monte Carlo tries to provide a more stable algorithm that allows you to find the constant in between. Geometric hierarchies are nearly optimal, but they are not the geometric hierarchies that you are used to work with. They come from this optimization. They have specific values that may depend on toll as well. So for ETOs, I included in the slide some assumptions and some results with respect to the densities that I don't have time to comment on, but I'm going to share the slides so if you have questions you can come to me. Again I emphasize here that these files actually depend on derivatives of the final solution with respect to initial conditions. So this is what really makes the method non-adapted to the natural filtration. So let's say this is multilevel in that context. This is what I meant by strong convergence and the accuracy result. This is what I meant by the work optimality result and I included a little bit more in terms of complexity. But OK, CLT2. So last thing before I run out of time, a little bit of an activity with respect to M and I just touch upon the classical expected value computation using just Monte Carlo just to close it up. So what's the point here? You want to compute the expected value of X using IID sample averages and you want to prescribe again this delta which is the confidence level and you want to prescribe also the toll. So what is the point here? The point is that usually you do this approximation based on the central limit theorem. That rule, that will be a rule of the following type. Essentially you are going to prescribe toll and you are going to prescribe delta and then a number of samples M will come out. It may not be so reliable when you have X that has low decay and density, it's essentially heavy tails. So for that reason we provided a rule that is a little bit more robust and essentially combines the central limit approximation with other type of bounds. Let me walk you through a little bit before I finish. I have 10 minutes, right? Yes, okay. So this is going to be a little bit intense then. So this is a general stopping rule, right? Based on the samples you are going to estimate the probability of failure and based on that we are going to accept or not. So the probability of failure actually is the probability of your approximation not being close enough to toll with the prescribed probability. So this event of failure having a value larger than delta essentially. So again as we described with the central limit theorem before the use of variance constraint you can actually, if you believe that your sample average, your sample variance is actually accurate enough you put a rule of the central limit type substituting sigma here which is a true standard deviation by its approximation and if this actually is accepted then you stop and you give your approximation. These may have issues, okay? So consider the following Pareto distribution. The density is decaying like X to the minus 4.1, right? For X larger than 1. And you see here the behavior of this quotient, this quotient is the probability of the actual failure versus delta the one that you wanted to impose. So this quotient is larger than 1, this is bad news, this quotient is less than 1 or equal than 1, it's okay. So you see here tolerance versus delta, right? So here it's like having more and more confidence in your computation and here is more and more accuracy as you go to the left, right? So as you let tolerance go to zero essentially for delta fixed things are going to be fine but in the regime where delta is say a few percent and sorry, when total is few percent and delta is sort of stringent and that's natural for engineering computations this is not going to be nice, it's going to be actually in the red and it's not going to give you the prescribed tolerance. So there are many results in the literature where you fix the delta and you let total go to zero. What I claim is that essentially those results overlook this regime which is important in the applications. Okay, so the whole point here is to penalize this fact because we really don't believe sometimes that we have such an accurate estimation of the variance to uncover this what we do essentially is to use very same type of bounds that are essentially sort of better control rates on the approximation in the central limit theorem and combining those with the sort of H-worth type of expansion which is something that you would use when you have full confidence that the tails are going fast essentially to zero then you get sort of an improved rule and when you go above this rule you actually see in this example that you can achieve the right prescribed failure rate with almost no change, almost no change you see in the computational cost and the cost here is just determinants of samples because we are working with just Monte Carlo in this context. And this can be done also for another kind of not so nicely behaving random variable and you see again that for small delta and relatively few percent total this is a bad case but you can fix it using this algorithm and without incurring too much in an increase of the work. And finally this is an example where we actually expect the usual approach of central limit theorem to work nicely and it does actually the work here is such that the failure probability is below the delta but this algorithm also achieved that with that an increase, a sensible increase in the work. So in the nice cases this behaves well, in the worst cases or the bad sort of cases this will have better control on what is going on. Okay so that's a summary and I think I get to my end that's remarkable actually. So thank you. There any questions or comments? Yes, Amar. So I would be interested to know if we could apply this criterion for example in the case where X is related to Rar event. So that is a Bernoulli random variable typically which is 1 and 0 but 1 with very small probability. What? That's not the case we looked into. I think that needs a measure change. I don't think one can apply it directly here. Just using some simple averaging of your output. So our robust is the criterion to... I mean what I'm saying is that it can be used but provided the first you do a measure change. If you don't do a measure change all the estimates are going to be wrong, right? You mean wrong inaccurate in fact? Absolutely. The whole thing is completely reliable, right? So Rar event asks for a measure change to start with and then if the measure change is sound, the important sampling is sound then you can use things like this. But I wouldn't recommend using this in a Rar event without looking into an important sampling algorithm first. It can be combined with the same... That's what I'm saying. But I wouldn't recommend it to use it without the important sampling step first. I have a second question perhaps. Because again that regime is different, right? I mean that regime is like for example 0, 1 with a very tiny probability of achieving 1. I mean here it's like unbounded, right? You can get unbounded variables with a certain probability and this is the case, not that fast. So you can achieve them with the probability that it's not so negligible and this is what gives you trouble in the estimation of the basis. It's not the same regime as the... Have you looked in the direction of concentration of measure inequalities like we know... Well, that's also not applicable here in the sense that for concentration of measures at least if you use classical large deviations you need exponential moments. These guys do not have any, right? I mean we use them in the Cherno for the Taulib but here you don't have an exponential moment. So if you have a concentration it's too little. A short comment. Maybe the applicability of your method is connected with a kind of ratio between the expectation and the variance. Because the case of Fourier events is that 1 over square root of p goes to infinity and p goes to zero, something like that. The relative value of the expectation... Yes, yes, but I mean that's also what defines the regime of red events, right? Yeah. But in this case again I emphasize that is the tail that is creating the issue and we create sort of something that works more or less uniformly with respect to that and does not address at all the question of the red event. The red event should be addressed in a complementary way. At least it's my humble opinion. Other questions, yes? Have you looked or is there any chance to extend this framework to the case of PySYS adaptivity? I'm thinking of SPDE where you would like to have a discretization that has some local adaptation in space. But this is not about this last part, it's the previous part. Yes. Well, I'm kidding. Okay, okay, okay. This is about the... Again, the adaptivity would respect to... Pathwise. Pathwise, okay. So what we discussed was pathwise but in time. Yeah, but the time step is fixed. No. Independent of the... No, no, no, no, no. The time step is delta T of T and omega. Delta T is fully stochastic. So it cannot have two features that come with omega, pathways. So yeah, that's actually the most general but you only use it if you're obliged to. That's my advice. That if you can use, for example, the terministic type of activity, that is cheaper. But you need some prior information about the problem that you're solving. Right? Okay, thank you. Did you have any photo questions? If not, so let's thank again the speaker. Thank you. Thank you. Thank you. Thank you. Thank you. Thank you.
 Thank you very much for coming with me today.  And I'll give two lectures.  The first one will be related mainly to multilayer Monte Carlo and the second one we will touch  upon multi in the expressions of the things that I'm going to say today, but not in an  adaptive fashion as we will see.  So first go then to discuss a little bit about the activity within multi-level Monte Carlo.  Let's see how this goes.  So the plan is as follows.  I will, well, these are some pictures from the place I come.  This is a chaos actually.  It's a King Abdullah University of Science and Technology.  Some of the students are here today.  And let's see.  This is a plan.  So most of you are very well acquainted with Monte Carlo and possibly multilayer Monte  Carlo as well.  So I will recap a little bit just to set up the notation at the beginning.  So we will speak the same language.  And then we go into different versions of our activity.  As you would see, there are many ideas for adaptive approximations in multilayer Monte  Carlo that we can see.  And for that reason also, I will touch upon different models.  I mean, one of them will be computations with random PDEs, that is PDEs that have random  coefficients.  And the second one will be diffusions.  These are etostochastic differential equations.  And the third one will be pure jam processes, say stochastic reaction networks, whose trajectories  are piecewise constant.  So I think this will motivate a little bit different ideas that we will have in this.  And finally, we will discuss a little bit of the contribution multilevel and optimal  hierarchies.  Hopefully, if the time allows, I will also discuss optimal stopping with non-asymptotic  rules within Monte Carlo.  So I'll try to really be adaptive in this process.  So let's begin with our first example.  This is a linear elliptic PDE that is very popular in the context of random PDEs as a  model case.  So here, the coefficients a and f are random processes, stochastic random fields.  And mainly, we'll concentrate on the effect of a, because a is with the two, right, the  one that acts non-linearly on the solution u.  And out of the solution u, we are interested in computing a statistic, which is real value.  Again, as an example, in this lecture, I'm going to go through the computation of the  expected value of a functional of the solution here depicted by psi of u.  Psi could be linear or non-linear.  We have some restrictions on smoothness on psi.  But for the moment, think of a linear functional that is the magnetic acting on u.  Again, this is not the only thing that you can compute with Monte Carlo, but this isn't  the one that we will be focusing mostly through the lectures.  All right.  So the goal is set.  We want to compute this expectation.  Now this expectation involves the solution of a differential equation.  Seldom, do we have the actual exact solution of u of our disposal?  So instead of u, we're going to approximate it by uh, say, where uh comes through a discretization  of this equation.  In this case, of course, it will be only space discretization.  But in general, you may have time as well or even other parameters.  So think that you have prescribed the use of your favorite method.  In this context, you could define it elements, find a difference, find a volume, so whatever.  And that discretization is indexed by the parameter h.  All right.  So h will govern in two aspects of the solution.  One will be the approximation characteristic of uh to u, but at the same time, the cost  of evaluating uh.  With the assumption, of course, of consistency, meaning that when h goes to zero, uh converges  in the proper space to u.  And at the same time, the cost of evaluating uh goes to infinity.  Right?  Otherwise, we will let h go to zero and get the exact solution.  Now, since we were interested in this function actually, which is a real value random variable,  we will approximate, I mean, the notation denote the approximation by gh, where now  h is the discretization parameter.  All right?  So that's hopefully a natural notation.  Let's go through the error then.  The difference between what we would like to compute, which is the expected value of  g minus what we actually compute, which is a sample average based on IID samples out  of the inexact quantity, which is gh and not g.  Right?  So we split upon two main sources.  Right?  One, again, the trick is to somehow subtract the expected value of this approximate quantity  to split the difference between this and this into two things.  Right?  One that we call the bias error, the other is a statistical error.  And for the bias error, we have essentially an estimate of this kind, which is coming  directly essentially from the deterministic analysis, right, point-wise, and then integrated  upon this expectation.  So essentially what we have is that if your method has order w here with respect to h,  which is the discretization parameter, when essentially the constant that comes from the  deterministic analysis will be random, of course, because it will be depending on a and  an f.  But once you take the expected value, you get the average of such a constant, which is  no longer stochastic, it's just a deterministic constant.  All right?  So this analysis holds as long as there are no changes of order with respect to the randomness  in the problem, which is a reasonable assumption.  But if that happens, actually, this will not be the correct way to draw things.  But again, this is an actually way to explain how this works.  So let's again consider that this w doesn't change with respect to omega, and therefore  we just get a constant here times h to the w.  Okay, so this is more or less a standard discretization error.  Second term is the standard statistical error, right?  And this error essentially is mainly governed by the number of samples in.  And I say mainly because the random variable that you're actually sampling depends on h,  right?  So this is an unbiased error, correct, because we take expected value here and then some  sample of the same random variable in an ID fashion.  So we only need to control its variance, right?  So and that is motivated by the central limit here, essentially properly scaled, right?  This will be approximated by a normal distributed random variable.  And this allows to write then later the constraint of having this error smaller than a certain  quantity in probability, you know, as a constraint in the variance, right?  So as we can see, again, we have two contributions here, one coming from h, the other coming  from m.  And as long as the problem is a rational well-posed problem and there's nothing strange, these  variants is uniformly bounded with respect to h, right?  And therefore they claim that m is essentially the quantity that governs the statistical error,  okay?  And if the problem has a singularity of some kind, maybe this is not completely true, but  I mean, again, in the examples that we're going to visit, this variance is going to  be bounded uniformly with respect to h, okay?  So two sources, one is governed entirely by h and converges algebraically with respect  to h.  The w here is determined by the numerical method that you use and the regularity of the problem  that you're solving.  And this m is very robustly coming into this analysis, right?  We only need to ensure that this variance is bounded, essentially.  So now that we have the two main errors identified, we try to see what is the corresponding complexity,  which is the computational work that is needed to achieve a certain accuracy with Monte Carlo  in this context for a prescribed confidence level, okay?  And to do that, essentially what you do is to minimize the work subject to the total  error constraint, okay?  So the total error constraint, again, improbability, because this was essentially an analyzing probability,  is the sum of the bias plus the statistical error.  The bias was governed by this, the statistical error governed by that.  And the total work, naturally, is what?  Is the product between the number of samples that you have to solve for, right?  Sometimes the expected cost per sample, okay?  And again, this is an example, right?  And the example that we have illustrated is a linear elliptic differential equation,  right?  And in that case, what we have to do is to solve a linear system every time we need  to create a new sample of the random variable.  So if v is the dimensionality of the space that you're solving, right?  I mean, could be 3D, for example, 2D, 1D.  Then h to the minus d, if h is a discretization parameter, is essentially the number of degrees  of freedom that you have in your discretization.  So h to the minus d is a degrees of freedom.  And then to the power gamma is essentially the cost that one incurs by solving the system,  okay?  So gamma equals one will be an optimal solver, and gamma equals three will be the Gaussian  elimination, okay?  So gamma can raise between one and three.  And of course, h to the minus d is the number of degrees of freedom in your system.  Is that clear?  Okay?  So this adapts to, again, to the type of solver.  If you use one type of solver or another, this gamma may change, okay?  Let's see what else is to be said.  Yes, okay.  So I mentioned the following.  I said that this cost was understood in the mean sense.  And why do I say that?  Well, if actually you have an iterative solver, right?  And the spectrum of the operator that you're solving with is random, right?  Because the a is the feasibility coefficient, was random itself, the number of iterations  may depend on the realization, right?  So what you're really interested in is the mean cost, not just the cost.  I mean, the mean cost will be representative.  And there will be variations with respect to each of the realizations.  Is that clear?  Okay?  Now, if you use a direct solver, this will not be the case.  But okay, I just want to emphasize a little bit on that.  Now, if you actually do the minimization of this bound and subject to this constraint,  you get this classical cost with respect to the tolerance that you're imposing.  Okay?  In the context of the central limit theorem, the confidence level is only acting on the  C2.  So it really plays a role in this constant.  I mean, in this exponent, it's just acting on the constant.  And that's why you don't see it.  But again, everything now depends on the accuracy you want to achieve and can be split into  two contributions.  One that is told to the minus 2 and the other is told to the minus d gamma divided by w.  Now the first part, told to the minus 2, is the classical cost in terms of number of  samples of a Monte Carlo method.  This is essentially what is needed to make this less and tall.  So that will be the cost that one will have to pay if each sample cost me over 1.  The second part of the cost is coming from the fact that I need to make h sufficiently  small so the bias constraint is satisfied, meaning that this term is actually less and  tall.  Okay?  So essentially what you see here is the number of samples times the cost per sample, which  is essentially what you were seeing here as well.  Okay?  Now, what Monte Carlo tries to do is to go beyond this bound.  So what we will improve somehow is this part and this will still be a barrier that we will  not be able to cross at least with Monte Carlo methods.  So just to give the notation again for multilevel Monte Carlo, let's recall this construction  of Henry and Giles.  Let's take a constant beta for subdivision.  For all purposes, think of beta equals 2 for the moment.  Okay?  So what we say is the following.  Let's just take in a single h that was meant to satisfy the bias constraint.  What we're going to do now is to take a number of such h's.  In principle, if you come from analysis, from numerical analysis, this is a natural way  to proceed because every time you are actually trying to solve a differential equation, you  try with a certain mesh.  If you're not happy, you subdivide it right and you test the stopping criterion.  And if you're not happy, you subdivide again until you meet the stopping criterion.  Normally, what you do is that once you get to that level that you accept the error to  be smaller than your prescribed accuracy, you throw all the other things, right?  You don't care, essentially, and you just work with the last accepted guy.  Here on the contrary, you're going to keep those for a certain purpose.  And the certain purpose is to produce an effective control variant.  So what one is doing here, again, let's go through the definitions first, is to introduce  these delta operators that are differences between consecutive approximations.  So essentially, this, remember, is a functional evaluated using a mesh of size hl.  And this is a functional using a mesh of size hl minus 1.  This is the difference between two random variables.  And we are assuming, of course, that they are not independent.  They are going to be sampled according to the same omega.  That's fundamental in this construction, because that is what will make this difference  actually go to zero, essentially, when l, bit l here goes to infinity.  That again comes from consistency.  In l2, this is what we're going to mainly need for this to be working fine.  So let's go again to the construction.  The construction of these hls is done via a geometric discretization, meaning that you  start with a given mesh.  If you're not happy, so divided by 2.  If you're not happy, so divided by 2 again, by 2 and so forth.  That creates this hl based on h0 divided by beta to the minus l with beta equals 2.  Is that clear?  So now that these deltas have been defined, then one rewrites.  Remember this was the guy that we were trying to use as an approximation to this one, so  that the bias between these two was less than tall.  Now this is somehow the random variable that one is going to apply control variance to,  and before doing so, one rewrites this expectation as this other sum.  So this is an exact representation, and the sum is running overall levels of approximation  from zero to the level l capital that satisfies the bias constraint.  Yes?  Now instead of using multicadload here, one instead uses multicadload independently for  each of the terms coming into this sum.  All right?  So, and that's sort of a scary proposition because the first time you see here, this  is saying at the beginning I only had one expectation to compute.  Now I have l capital plus one.  How is it going to work?  Well, okay, before we look into that, let's see how this looks.  And as you see, this, which is actually the estimator for multicadload, looks precisely  like the sum of l capital plus one multicadload estimators, where each of them now is sampling  something different, each of them is sampling these deltas.  And of course, remember that for l equals zero, actually, this is not the difference,  it's just g evaluated on the course's approximation.  Okay, so you can think of this as, you know, g is h0, right, if you think of hierarchy  approximation, and then the other things are just corrections, right, that are being added  to the course's level.  Now that we have this, how do we motivate, at least first intuitively, why is this going  to work better than the original multicadload approach?  Well, let's see, the bias is not going to be changed.  We said, essentially, if this L is chosen so that, you know, this bias is less than  told, what we are rewriting, essentially, is this quantity.  So if you're going to actually improve this sampling is because we're going to pay less  in the samples, essentially.  The difference will not come from the bias.  Yes, please.  So how do we choose for each level?  Oh, this is coming.  This is coming.  This is coming, don't worry.  This is coming.  Okay.  Now, the variance here is the variance of multicadload, right?  So what did we say that the constraint, you know, the statistical constraint looked like?  Well, essentially, thanks to the construction, based on independent samples, the variance  of this estimator is just one over ML capital of the variance of this GL capital.  Now, you think that the GL capital is more or less close to G, so more or less you can  expect that this looks like the variance of G divided by ML, right?  So this is essentially saying that in the deepest level of approximation, the capital,  the one that satisfies the bias constraint, you need to use total to minus two samples,  right?  And that is what you expect from multicadload.  If you, on the contrary, use multi-level multicadload, right?  The variance for the estimator, again, by the construction of independence, not only  between the samples that we use on each level, but among levels, right?  We get some of variances across levels.  The first variance, the one that is coming from the level zero, actually, looks very  much like this term, right?  This is not a difference.  This is not a difference as well.  But the other terms are actually one over M, right, times variance of differences.  And what's the point here?  Remember, as little L gets larger and larger, we hope to have convergence, okay, just by  consistency.  Now, in which sense we expect to have convergence in this setting, what is most useful is to  consider not only you have convergence in the mean, but also in L2, which means that  these variances of differences are going to zero, right?  This is where multi-level multicadload becomes effective.  And if they actually become very, very small, it means that for really deep levels when L  is large, right, each sample cost me a lot, but the variance is correspondingly small.  So the number M that we have to use is really small.  And this is where you get real advantages, because you are sampling much less in the  most costly levels, all right?  Now, you have to check some conditions usually to make sure that the choice of the level  zero say sound, but again, the whole thing can be put in terms of these inequalities  that again indicate that L2 is very relevant in the choice of the multi-level parameters.  So far so good?  Any questions?  You said there was to be correlation between the different levels.  No, there's not.  Actually, no, they are constructed totally independent.  Okay?  That's important.  Yes.  From the independence in the construction, you get immediately this sound.  Okay?  Otherwise, you will have extra terms, all right?  It's not that it's impossible to do that, but the construction here is done on independence  in some chance.  Okay.  So cartoon type of approximation essentially, we are somehow moving samples from sort of  deep refined levels that are costly into levels that are cheaper and less refined.  Okay?  Okay.  So if you want to go beyond that intuition, you need some assumptions, right?  And these assumptions may or may not hold in your problem.  You have to verify that they are satisfied and sometimes you may have precisely this kind  of convergence and some other type, I mean, depending on what problems you have.  So here we are again going to concentrate on this algebraic type of convergence with  respect to H. Remember, H was the parameter of discretization.  You can think of delta x, delta z, delta y in your linear elliptic PDE.  And we were assuming that the HL was essentially H0 times beta to the minus L, right?  A geometric subdivision.  So this first equation here is nothing different from what I wrote before, right?  It means that the error in the mean is H to the w, essentially.  You see?  So nothing is different.  It's just a rewriting of what we have before.  Now again, this w, which is the exponent of the weak convergence, depends on the assumptions  on your problem and the numerical discretization that you're using and has to be found for  each problem, right?  The variance assumption again is indicating some kind of error rate, but this error now  is measured in L2.  And instead of taking a difference between g and gl, now we are taking, remember, the  difference between gl and gl minus 1 when we put it inside of the difference.  And we have a strong rate of convergence that is s, which is not going to be w, right?  In general.  So for the last one, we have the cost, the mean cost of evaluating a difference on the  level L, right?  So this is delta gl.  And this again entails the cost of evaluating gl and gl minus 1 because we're going to take  the difference.  Usually when assumptions that gl is the one that dominates, so this goes like that.  And this is nothing from what we said before again, because before we said that the cost  was H to the minus d gamma.  And again, this is the same thing.  Just rewritten in terms of this beta and this L. I emphasize again that the constant  c, w, and s also depend on the problem and have to be found.  And remember, gamma may depend also on the solver that you use, not just on the problem  and the discretization.  So if you change the solver, gamma will change even though you keep the rest of the things  constant.  As an example, for a smooth linear PDE of the type, linear type that we discussed today,  if you use multi-linear piecewise continuous finite elements, we will have 2w equals s  equals 4 and of course, depending on the solver, use gamma ranging from 1 to 3.  So what's the total work?  And here we start getting close to your question about the m's.  The total work will be found in the same way as we did with Monte Carlo.  We have L capital plus 1 levels.  On each level, we have to use ML samples and the cost per sample.  It's the work of delta GL in a expected sense.  So what do we do?  Well, remember, we have to find ML.  So if you fix L, essentially, you can write an optimization problem and I'm going to go  into that a little bit more when I go into multi-index Monte Carlo, but we can just say  it again.  I mean, remember, if you fix L, the bias constraint is not going to be important and  what you want to minimize is minimize with respect to the ML's from L equals 0 to L  capital.  What?  The total work.  And this is the sum from L equals 0.  I hope it is possible to see this.  There's no reflection or anything like that.  So L capital ML and this work delta GL.  And then, of course, there's a constraint.  The constraint is the statistical constraint.  The statistical constraint, again, is the variance being less than the total square.  So what is the variance again is just 1 over M0 say variance delta, actually G0.  Just write it differently.  Plus the sum from L equals 1 to capital of the variance delta GL.  And this is less than total square.  Of course, there's a constant here in between, but I'm not writing it.  The constant comes from the confidence level that you want to impose.  Okay.  So the control variables here are ML, capital plus 1, right?  And here you see how they act.  So you can just write the last of relaxation and you realize that you can find these ML's.  And ML's actually are going to be, at least if you make a real value relaxation first,  this will be total to the minus 2.  And then there's a constant here that is coming out of the levels.  And this essentially is the variance of delta GL, I think.  You may correct me if I'm doing this wrong.  And then this is the cost, essentially.  And this is for root here, if I'm not mistaken.  Okay.  So you can first write this lagrange of relaxation, solve it, and then say round it above, I mean,  take the ceiling of this number to make sure that everything goes well.  Okay.  And from there, you get to that formula, essentially, assuming that the cost is not trivial, essentially,  that you're not getting less than one realization per level.  But that is a technical assumption.  You can stay with this for the moment.  So what do you see from here?  You see that after you substitute the optimal number of samples, the cost only depends on  L capital.  Because these are problem, say, and discretization-based constants.  So once you give the discretization method and the problem, of course, this work L and  the variance L, essentially, are given.  Maybe you don't know them, but precisely a priority, but they're given.  And this is a sum of positive terms.  So first thing that you observe is that when the tolerance goes to zero, L has to go to  infinity.  Right?  Is it clear why?  You have to satisfy the bias constraint, right?  So the bias constraint, essentially, was h to the w less than toll.  And L is coupled directly with this h to the w.  So L, essentially, goes to infinity like a logarithm of toll, absolute value, say, that's  a ballpark.  You may have log log terms on top of that, but really it's log toll.  So this is going to infinity as toll goes to zero.  So this sum can only do two things, either stay bounded as toll goes to zero or diverge.  If it is diverging, it's because this work is somehow growing faster than what the variance  is being able to converge, right?  But if it is actually finite, which is essentially the case where this variance is dominating  the work somehow, we have a very interesting case.  But it's sort of not intuitive the first time you see it as well.  It's the fact that you get a cost that is essentially all of toll to the minus two.  And if you remember Monte Carlo, when we discussed it first, toll to the minus two was the optimal  cost with Monte Carlo when each of the samples, essentially, in the mean cost me all one.  So essentially, this means that you don't feel the cost of discretization even though  you're refining the levels, right, infinitely much as toll goes to zero.  Is that okay?  In principle, not, right?  So somehow it means that even though we still have to sample of those at the level L capital,  we only sample so few, so little samples so that that cost really does not influence the  rest of the cost that come from the other levels.  Because we still have to sample in that way.  That is clear, right?  Okay, so now that you see what is coming, let's see the result with these assumptions  in mind.  L capital, as I told you, is determined by the bias constraint and essentially is linked  to toll by a log factor, right, where the log essentially comes connected with a weak  rate of convergence and the logarithm of the subdivision multiplicative constant beta.  And the samples ML are chosen optimally like we discussed before just by solving that Lagrangian  relaxation problem.  And after substitution, essentially we get this fantastic case, toll to the minus two,  where S is larger than the gamma, right?  So essentially both the work and the variance are behaving with respect to L in an exponential  fashion, right?  So the work is growing with the D gamma rate and the variance is decreasing with the S  rate.  So as long as S is larger than the gamma, this sum is convergent.  And this is what explains the first result.  The second is a limit case where essentially they outmaneuver each other.  So essentially they balance each other.  You get something that is constant here.  So the total result here is L on this sum, L times a constant, and then you get L squared  and L was log toll.  So that explains the second line here.  And in the last one, well, the work is growing faster than the variance is converging and  instead of getting something that is convergent, you get some algebraic deterioration.  But still this algebraic deterioration is better than what you would get with Monte Carlo.  So you're decreasing from D gamma here in the explosion to D gamma minus S.  OK?  So you're still exploiting a little bit of the strong convergence.  I mean, the convergence in it too that one assumed.  So OK?  So far for the introduction with the uniform methods.  Now the question that we want to more or less touch upon is how to extend this Monte  Carlo into non-uniform adaptive discretization settings and whether it's worth it.  OK?  And there are several ideas to introduce non-uniform discretizations or stochastic levels, if you  wish.  So let's begin with the...  OK, well, this is...  I put just to include in this slides and I already gave up trying because, I mean, the  pages are not enough actually to describe so many things.  And I rather just direct people to look into mics.  Himongu said to keep up with the growth of the community.  And this is actually what you get the best references in the...  List of references in the area, so I'm not going to talk too much into that.  So let's go to the activity then.  And in this context, again, we're trying to approximate the expected value of an output  quantity.  In this case, I'm switching the example.  I'm going to stochastic differential equations, itto differential equations.  So in this setting, each of the paths evolved in finite dimensions with respect to time,  right?  So you can do this with the sum of the paths to a final value and then in that final value,  you may evaluate the function g here that comes in this example, like a put option in  finance and then you have to integrate this g versus the PDF of the final time.  OK, so you can do this with Monte Carlo just by sampling many, many paths and again averaging  as we discussed before.  So here are some references again.  I'm not going to stop and I'm going directly into the problem formulation.  Here is an SE.  This SE, if you have not seen it before, it's not a problem because you can think of this  discretization via Euler-Marujama that is coming just in a second.  So x is something that is finite dimensional now, it's not a function as in the first example.  It evolves in RD.  And g is a given function so that our goal is to compute the expected value of g of xt  and so the evaluation at the final time of the solution of an SE.  Again we are given an accuracy constraint tall and a prescribed confidence level and  everything is driven by k-dimensional venerable process.  All right, this has plenty of applications that I will not stop upon and what we're going  to do is through this adaptive Monte Carlo and actually see that there are several versions  now of adaptivity.  The first one we did was with respect to the weak error but later we did something that  controls the strong error and yeah, let's go into that later.  So first discretization algorithm.  This is just Euler-Marujama so if you have not seen SEs before, this is a new term that  comes in.  So remember it's coming through these.  If you just think of the first part, it's just forward Euler right based on this function,  the drift function a and the second is just the same type of discretization.  You see you have to evaluate b on the left point, essentially the current point x bar  but you multiply times these random variables.  These random variables are normal distributed with mean zero and the variance is just the  delta t that you're taking.  So the delta t does not need to be uniform at all but yeah and again weekly you can just  scale them to normal zero one random variables and the multicolline estimator looks exactly  as we discussed before.  It's an average of IID samples but the IID samples are not from the exact process but  from the solution of coming out from this forward Euler.  Forward Euler of course is not the only discretization that one can do in this context but it's the  one that I'm going to exemplify.  Okay so as we discussed before this can be split into tolerance for time discretization  and tolerance for statistical error and this is a particular case as we discussed before  the cost per sample here is essentially proportional to the number of time steps.  The number of time steps essentially is proportional to one over a ton.  This coming from the weak convergence rate one of Euler which is the same as Euler-Maluchin  in this setting and this brings us to the complexity told to the minus three for Monte  Carlo.  Okay okay and now you want to introduce adaptivity.  So what is adaptivity in this context and this context adaptivity means the ability  to choose the time steps non-uniformly and this non-uniformly can depend also on omega.  So the time steps such as t0 is 0 and the final time step is t capital are not going  to be equally spaces and they may adapt to omega because the dynamics of the solution  may depend on omega in a dramatic way maybe.  So in this context then we will create realizations that depend on x bar, I mean here denoted  by x bar and depend on this mentions delta t of omega.  So again why do we do this?  Well sometimes there are non-smoothness on A or B that ask to, I mean to deteriorate  the usual behavior of the uniform method so ask essentially for the adaptivity in a natural  way and what do we need to be able to produce such delta t?  Well you need, at least with our approach, you need to rewrite the weak error in this  other way essentially, you need to introduce an error density which is a function rho w  that does essentially not depend on delta t so that the error, the weak error can be  split essentially in contributions that are coming out of each of the different discretization  steps.  You see here if you discretize over the time steps, this is a sum over time steps of essentially  delta t n square times rho w.  Now the fact that delta t here is inside of the expectation is allowing for the choice  of stochastic time steps because rho w may still vary dramatically with respect to omega,  that's the point.  If this rho w will not be varying dramatically with respect to omega then it's cheaper actually  to take a non-uniform but deterministic discretization with respect to delta t, is that clear?  There's another head to do this fully stochastic so whenever possible you avoid that.  Now the funny thing is that in this context there's not only a weak density but there's  also a strong density.  You can actually make an approximation and you see here the power is one but here the  power is two, you can actually write the L2 error here in the same fashion for Euler  where now this rho s is not exactly this rho w, this is still something that you can approximate  with computable quantities and it's able to distribute the strong error across the different  intervals that come into discretization.  In principle this is not obvious because when you square this you may get a lot of cross  terms.  So to be able to tell that actually the diagonal terms are the one dominant this needs some  understanding.  Okay, but what is the whole point now is that once you actually produce these error densities  you can produce single-level adaptive algorithms.  So single-level adaptive algorithms will try to find these time steps so that say for example  a bias error is satisfied or even say this type of L2 error is satisfied.  And what does it allow then?  It allows for the following construction.  If you want to create a multi-level adaptive method what you derive essentially the discretization  with is with a sequence of tolerances.  So if you have a reliable single level adaptive method the only step that you need to go from  there to a multi-level adaptive method is just to create a sequence of tolerances and  to say that the first level will be driven by what will be defined by discretization  of my single level adaptive method so that they reach that accuracy.  The next level is again refinements on those that achieve the next accuracy constraint and  so forth for all the levels.  So you see now each of the levels does not have a prescribed H. This notion is gone.  Each of the realizations may use a completely different mesh that adapts to the characteristics  of the omega that you have in mind.  But still one will be able to extract the good properties of the variance convergence,  the mean convergence and the rates of the explosion of the work with respect to the  level.  Okay?  So what is optimal at least to show results on the algorithm in this context is to refine  according and to stop according to these simple rules.  Essentially what you do is to compute solutions and error indicators are n for each time step.  And the error indicators are n in the context of this error density defined essentially  as rho n times delta tn squared.  I mean and rho n could be the weak error if you're controlling the weak error of the strong,  if you're using the strong approximation that's not important.  But again now this is defining the contribution from times the n in a given realization to  the total work.  And what we do to stop the discretization is to verify whether the maximum of these guys  is or not above this.  And this you see is the tolerance that we prescribe for the bias constraint which is  what we want to achieve divided by the expected value of a number of time steps.  So this essentially is the mean error that we would like to achieve on each interval.  That turns to be optimal essentially to equidistribute the error along intervals.  One caveat though, this n is not known a priori, right?  Because this is the number of steps that you need to solve.  So these have to be found on the fly.  But this is a constant and this is another constant that are essentially probably independent.  So what you have to remember here is that this is going to be not the same as this C  of refinement.  The refinement criterion is a little bit more stringent than the acceptance.  So essentially if you see that you're not going to stop because some of them are going  to be larger than that, well then you are going to refine many more than just those  that not satisfy this constraint.  Now another technical point in this construction is the following, if you have a refinement  and has not been accepted and you mark some intervals in time for refinement, you need  to produce essentially the corresponding interpolation values for the delta W's that drive your discretization.  And these, you know, not to introduce any biases, have to be consistent with the points  that you already sampled.  And this is done according to Branian Bridges, right?  That is the correct way to do this discretization.  So this plot essentially shows that, I mean, that for a given level, for example, you use  the blue discretization, right?  And this entails evaluating W at all those points.  But now you need to evaluate W in points in between.  And of course, they are not going to be linear interpolants of the previous values, but they  will differ upon.  And this difference is actually given by the Branian Bridges, okay?  So it's just conditional Gaussian sampling.  That's nothing more.  Okay, any questions?  No?  Okay, so what was the main idea again?  Let the tolerance in the delta W algorithm define the hierarchy.  That's essentially what you have to remember.  And this hopes, this works in any problem.  If you have a sound error density of PEs or EZEs, you can just use this approach and  it will work.  Like I said before, adaptivity is useful when you have issues with A or B or the process  itself may have issues with respect to, say, boundaries where it is being, for example,  stopped or reflected that those actually ask for adaptivity.  You may have stiffness as well, and this is also nicely treated by adaptive algorithms.  So again, if you have a stiff drift, you may choose at least delta T.  If you don't want to go through this posterior global error control, if you just want to  drive things by stability, you just may take delta T to make sure that you remain within  the stability region of the method.  And that again, if the stiffness is just coming from A, it can be done quite easily.  So what is interesting here?  The point is if only when you do a posterior error estimate, you show things about the  discretization being efficient or not in the sense that the estimate for the error essentially  gets close to the error you want to impose.  But in this context of error densities, you can show actually, you show much more, you  can show actually that the algorithm stops, that asymptotically you have normality, asymptotically  you have accuracy.  And again, what is asymptotic in this case is when tolerance goes to zero.  There is no notion of H equal to zero because H is gone.  So the only parameter that drives the discretization here is the accuracy constraint, and the accuracy  is given by toll.  So when you look at asymptotics in this context, it's when toll goes to zero.  And of course, we have kept the confidence level constant in this discussion.  We are not letting the confidence goes to one.  So asymptotic normality, you also have asymptotic accuracy.  I will just tell you later what this means.  But it really means that when tolerance goes to zero, essentially the error is bounded by  toll.  And the complexity essentially is up to essentially a constant, you get the optimal complexity.  So for the ethos to classic differential equations, this one discussion that I have  not made, and I just put the warning here, it's actually working fine in the theory,  but it's not obvious.  The discretizations that you get are not adapted, I mean, they are adaptive, but they are not  adapted.  Adapted means that they are not adapted to the natural filtration of the problem, which  is given by the linear process.  Why?  It's because the density, and I have not told you maybe, no, these densities essentially  involve variations of the final value with respect to the current value.  So essentially, they may need, for example, derivatives of respect to the position Xn  of g of xt capital bar.  This quantity is not adapted.  It's looking to the future from, say, the current point n to final time.  So in principle, you can discuss whether these things are going to converge to the correct  solution.  They do.  But I mean, this is not true for any adaptive method that you can come up with.  And no one has to be worried, I mean, at least careful about it.  Okay, so that's why.  So if you do examples here, for example, you have a drift singularity, this is just going  to be simple in the case, in the sense that the stochasticity is really not relevant.  I mean, here we have, say, deterministic time for the singularity, and the mesh is going  to refine automatically around the value of alpha.  We are not telling the algorithm that alpha is there, but it does find it.  You see the mesh here is refined from the level, essentially, 10 to the 0 to 10 to the  minus 10, just to control the error.  Okay.  But again, like I told you before, this is sort of the deterministic type of singularity.  So in this case, it doesn't pay too much to let delta w to depend on omega.  Delta w is not uniform, but it really doesn't depend on omega.  So you can split the construction here in a first phase, where you actually find all  these meshes that you're going to use in all the levels.  And then you just run in the second phase, material multicolor increasing the number  of samples that you need on each level until you hit the variance constraint.  That is effective.  Okay.  And actually, what you can show is that the cost of the first phase to find actually the  optimal discretization in Monte Carlo is much cheaper than the actual cost to solve  the problem.  So that's again, actually, when you have deterministic type of singularities that you can find the  mesh very quickly.  Okay.  So in this case, you see an advantage from uniform to adaptive.  You don't see much more of a difference between uniform Monte Carlo and adaptive Monte Carlo,  because essentially what happens is that this singularity is not strong enough to modify  the strong rate.  It does modify the weak rate.  But if you remember what I said before, the complexity in Monte Carlo is mainly governed  by the strong rate.  I mean, there's this independence of the weak, of course, but in the inequalities that I  showed you before with these conversions of the sum, the strong rate was the one that  was opposing the weak, the work growth, not the weak.  Okay.  So, but if you actually go to stronger singularities, so now we let P actually be beyond one half  and, of course, you have to make sense that this is actually has a unique solution, blah,  blah, blah, but it does actually, as long as P, I think is less than one strictly.  So you can see here for P equals 0.5, that was the previous example, you don't see much  of an advantage between multi-level Monte Carlo uniform and adaptive, but as you start  raising the P, you clearly see the advantage between adaptivity and non-adaptivity.  And again, I have to tell you one more thing, this has been done using the strong density,  not the weak density.  So this is the latest work that we have for adaptivity.  It's actually on the strong, I mean, you use the strong density to rewrite this L2 error  and you now create the levels according to this strong density.  That's how this plot is being created.  One more thing, okay, the alpha here was uniformly distributed between one fourth and three fourths,  so actually the measures have to be stochastic, because the position of the singularities  is being just jumping around.  Okay, so I told you essentially how to do adaptive multi-level Monte Carlo approximation,  yes please.  Sorry, alpha is independent of w, that's right.  So essentially you set alpha, you run the things, alpha is independent of w, correct.  Yes, I mean, this was just a simple example.  So people can actually reproduce and see what is going on.  But yes, I mean, things can be more complicated because the singularity, as you're saying,  can be coupled with the dynamics, okay.  So it could be that actually the singularities appear because of some drift problem, right,  with the stiffness, or it could be that x is getting close to a boundary and the boundary  has an issue.  So you can think of the boundary being penalized and penalization is drift stiffness, right,  so yes.  But that will not be able to be written in one line essentially, so yes, at least this  motivates what you can gain.  Thank you.  Okay, any more questions at this point?  I'd like to then summarize what we did on this topic.  So I hope the idea of adaptivity in this context that I gave you is quite general.  You can use it for any problems you want essentially, as long as you have a sound single level analysis,  if I may say, okay.  The message is that it is better in the context of multilevel to use an analysis that will  go towards the L2 error, if possible, right, drive the analysis of the construction of  the levels according to that.  Of course, you need the bias error to control the number of levels you're going to use,  but at least the distribution of the shape of the meshes should be done according to  the strong error.  And actually that can be proved, actually, that is optimal to do so.  You just write a problem that is related to these, but instead of writing this part here,  you just write the error density times the, I mean, the integers of the error density  times the delta t's, you just write the whole Lagrangian relaxation, which now will involve  finding functions, functions delta tL, say, from L equals 0 to L capital, and then you  will see that delta tL is directly expressed in terms of the row strong, not the row weak.  The only, the last level is connected to the row weak, but all the rest are driven by the  row strong, and therefore that motivates directly to the use of the strong for the activity  that we found later, sorry, but, I mean, the first works we did was with the weak activity.  So how do we do with the time?  We split in between or?  As you like.  Okay, I can just draw the conclusions and it would take five minutes.  Let me put the conclusions and then, okay.  So as conclusions, I told you again that this way of thinking for multi-hebel and adaptivity  is quite general, you can use it anywhere you want.  You have to be careful, though, that when you replace the statistical constraint, which  means that the statistical error has to be less than tall with a prescribed probability,  computationally it's much more effective to write that into a variance constraint, and  this can be done as long as you have a central limit type of approximation.  The central limit does not apply directly in the multi-hebel-vontical setting because  you have some of independent but not identically distributed random variables, so you have to  use resort to some linear or failure type of result, and you have to verify again these  four each of the problems that you work with.  So you have to be careful with that.  Okay, so let's leave the puzzle random measures for the next one.  Thank you.  Five minutes?  Yes, sir.  Five minutes.  Thank you.  So we are going to discuss another idea for adaptivity in this second context.  The second context involves paths that are also random and depend on time, but as opposed  to the stochastic differential equations that we were given before that were etotype, these  paths actually can be solved exactly.  So these paths are given by piecewise constant functions, and the only random thing about  them is essentially the time where they jump, and which kind of jumps we have them.  These processes are very popular, they are macro processes of course, and they are used  in Xenium applications, but let's go a little bit more into them.  This is a logarithmic scale, this is time, and here in blue, red, and green are three  components of this vector that is changing according to time.  The numbers that are shown here involving are the number of species for each of these  three components, so you can think of three populations randomly changing inside of a  system according to maybe biological reactions or chemistry reactions for instance.  Even though these green paths do not show it, they are as discontinuous as the blue  ones essentially.  It's just that you have logarithmic scale and you may be fooled by the logarithmic scale,  that's it.  So that's a motivational example, this is a gene transcription and translation you have  here in this system, the component G that may react given G itself an M, M that actually  may react again and give M and P, and the P that actually goes into D when it meets  itself.  So you see that there are sort of reaction channels, essentially when this thing happens,  you know deterministically what is going to happen.  The only problem is that the time where this occurs is not going to be prescribed, a theory  is going to be random.  And let's think of a simple example, let's think of water say in a container and water  can be created from hydrogen and oxygen and it could also dissociate into hydrogen and  oxygen, so this is a reversible equation.  But the way the speed of these reactions actually happen depends of course on the temperature,  but if you fix the temperature out of this discussion, they depend for example on how  much hydrogen and oxygen is there because if there is no hydrogen, oxygen will never  meet hydrogen to create water, right?  So the likelihood of these events to happen depends on how dense these components are  in the system.  So and this can be mathematically described.  So what we want to do again will be to compute some given output quantity of interest at  a given final time up to a certain tolerance again with high probability.  So how would describe the evolution of such a system?  I told you before this was an X that lives on a lattice essentially because it's counting.  So the components of X essentially describe different species in the system.  So we count how much of the species one we have, how much of the species two we have  and so forth, okay?  So X is taking values on these lattice, right?  And it was according to time and of course omega.  Now we assume that we only have finitely many possible reaction channels like in the case  of water for example, we may have three components, one is how much hydrogen, how much oxygen,  how much water, right?  And we can only have water dissociating into hydrogen and oxygen and the other two meeting  to get water.  So when this oxygen and hydrogen meet, right, we assume that they can only create water.  So there's nothing stochastic about the fact when we know that they meet.  The question is when they meet and this is the stochastic part.  Okay, so when this thing happens, right, the J reaction type of in the system happens,  the only thing that we have is a change from X that was the state of the system into X  plus nu J.  In the case of this water for example, when oxygen and hydrogen meet to create water,  you have a plus one in the water and the corresponding minuses in the oxygen and the hydrogen.  Now, when does this happen?  Is given by the so-called intensity functions for the jumps or in this context also called  propensity functions.  These are non-negative functions, A, J, that are evaluated on the state X and give values  on R plus.  So that the probability of having a reaction of type J in an infinitesimal interval that  goes from T to T plus DT, given that you have the state X is essentially A, J of X times  DT.  So for instance, if you have no hydrogen in the system, you know that new water is not  going to be created.  So this A, J corresponding to that, that reaction will be zero in that interval until that situation  changes.  So again, we want to compute the expected value of a given quantity for the given of  several G and the applications are everywhere.  So I'm not going to stop on that.  So we have plenty of references and as well this is not complete, but we just wanted to  show a little bit the efforts on this multilevel version.  So what is interesting about these systems?  First thing as opposed to the previous discussion that we have is that exact algorithms are  available.  What do you mean by exact algorithms in this context?  You remember when we sample a la Monte Carlo to approximate quantities of expected type,  we always have statistical error unless the problem is trivial and has no variance and  then we are not discussing anything stochastic.  So exact in this context means pathways exact, essentially that with finite work you can simulate  a process that has the correct law.  So there is no bias error.  You can do unbiased simulation with all one cost, which means that already you have the  beautiful regime that we discussed before that was total to minus 2.  The cost is total to minus 2, the number of samples in Monte Carlo and then all one to  solve for each path in the mean.  So the whole discussion here is about constants.  So you do variance reduction in order to reduce the constant in front of the total to minus  2.  So that's the whole point in this context.  As opposed to the other case where we actually can push the complexity from total to minus  2 log square times total to minus 2 into total to minus 2.  Here it goes from total to minus 2 into total to minus 2.  And still you can get implicit gain.  So constants matter, that's what I'm trying to tell you.  So then if you have an exact algorithm, why are you going to introduce an inexact one?  Well, because the two things, because the statistical error is already there.  So having a bias that is zero in the presence of a statistical error may not be optimal  because the cost of simulating such exact path may be a lot.  So introducing a path that is a little bit cheaper and has some bias and the bias is  smaller than the statistical error will not change things dramatically from the error  point of view, but it will do it from the work.  So it is actually a good idea to actually introduce bias in a bias simulation.  Right?  Good.  So what is the problem then here is that essentially the cost comes from all these events that  are happening in the system.  So if you actually have a lot of molecules in a big reactor, you can think that there  are a lot of reactions happening.  And if you have to track all of them to do this exact simulation, the cost is going to  be horrible.  Okay?  So, well, essentially these interarrival times between jumps are exponentially distributed  according to, with intensity, there is a sum of all the AJs in the system.  And this can be a large number, depending on the system.  Okay?  So as opposed to that, then one would like to do approximate time stepping.  Time stepping that is cheaper but introduces a bias.  This has been done in, I mean, pretty much from the tau-lib, there is a forward-oiler  method to more refined versions of it.  And you have drawbacks, of course.  First is you're going to introduce a time discretization error that you didn't have.  Right?  And two, you may have non-intuitive, at least qualitatively wrong type of behavior.  Essentially, if you go back to the discussion that we had at the beginning, our process  is actually counting.  Counting molecules of different types in a system.  These are not going to be negative, right?  I mean, if you're doing the proper thing, the true trajectories will never cross zero.  But if you do these sort of approximate algorithms with a non-zero probability, you may cross  these things.  And you may want to avoid these events because they have an issue on the accuracy.  The global accuracy that you can get.  So there are ways to deal with these.  For example, the pre-lib is just to adjust adaptively the time step, you know, the control,  the one-step exit probability.  But we would like to always relate things with respect to the goal, which is to relate  the old type of errors that we have in this approximation with respect to the goal, which  is the expected value of g of xt capital.  So what we did in this line of adaptivity was to first create a single-level hybrid  algorithm that at each time step, adaptively switches between the exact path simulation  and this tau-lib sort of oiler in order to do what?  To minimize the computational work.  So why is this somehow needed?  When you're close to this zero level that you don't want to cross somehow, the tau-lib,  which is the approximate forward-lowered type of discretization, does have a cost that  is not being reduced somehow.  So you're not gaining too much from using it.  But the exact simulation essentially has a cost that I told you was essentially going  up with the number of particles in the system.  If you're going to get closer and closer to zero, that cost actually gets smaller.  So close to boundaries in this respect, the exact simulation is advantageous because it  doesn't really use bias and it's cost to do because the inter-rival times are relatively  large at least.  So this can be quantified and a rule that switches between these two things can be done.  So again, this is adaptivity, but now not adaptivity just with respect to time, but  between methods.  So you're combining two methods.  One that has sort of infinite accuracy or infinite order with respect to time discretization  and the other that has order one, which is essentially coming from Euler.  So that was one thing.  The other thing is that this idea can be generalized into multilevel Monte Carlo algorithms for  our control again.  And not only that, it's like you can do this according to the different channels.  So I told you before that when you introduce this, there are different J channels.  And with respect to these Js, what matters really is the propensity.  If the propensity is sort of large, you're taking a lot of events over a unit time step.  If the propensity is slow, these things fire not often.  So if you can actually split among the reactions, do a time splitting method and essentially  treat some of the reactions with the approximate method and some of the reactions with the  exact method.  This can be done as well, even in the context of multilevel.  And again, all this, what is going on here is that you're introducing biases and approximations  that have much, much less cost when compared to the unbiased algorithm.  This is why this gets useful.  So how do we go about it?  To motivate the different methods, the whole thing can be written through this course random  time representation as accounting algorithm.  Essentially the value, the number of species that you will have at time t, essentially  is the number of species you have at time zero, right?  Changed by what?  Changed by the number of times the Js reaction occurred.  And the number of times the Js reaction occurred is described here by a Poisson process whose  intensity is non-trivial, it's just the integral from zero to t of A J, the propensity of the  reaction J.  This is again a non-negative function and this is just a random time change on a given  Poisson process.  So this is a natural number that goes from zero to whatever and this is just the vector  that introduces the changes according to that reaction.  Is that clear?  So it's just bookkeeping, you see?  And if this A would not depend on X, this would be trivial to simulate, right?  You will just change accordingly the time in a unit of Poisson process and know exactly  when the jumps are going to occur in your system and that's it.  So to simplify this, when you want to do approximations, you just do a donor method, essentially you  freeze the value of A at the left point and approximate this increment by a Poisson random  variable with intensity A J times delta t, if you wish.  And the first time I wrote it, they put tau here and this is why they call it the tau-leap  method.  Okay, so first problem that you see here from this approximation point of view is to find  tau so that you're not going to cross these important boundaries with a given probability.  Delta is given.  So this tau will depend on X and delta and this will prevent you from jumping according  to this inequality.  What we did actually is to do some large deviation, turnoff type of bound to actually come with  a computable formula that one can use to ensure this inequality.  And then from this delta t, we can actually discuss whether it is worth taking the approximate  step or the exact one.  So and then one step farther is to go into a global probability of going out because  remember this discussion that I just made was based on one step, essentially.  If you're here, how you have to take delta t, so the probability of exiting in this time  is below a certain value.  Now what you really want is the probability of the whole path going out or not going out.  So you have to link the one step probability to the whole path probabilities and then link  these quantities into the effect of the global error that we really want to compute and approximate.  So this is the first step that is good because I mean in principle the probability of going  out for a path, if things are going to be done correctly, it's a rare event.  So you want to, instead of sampling that rare event, you want to compute quantities that  are sort of easy to deal with.  And this approach does it actually.  Okay, and like I told you before, then one can write a multilevel Monte Carlo based on  these type of paths, right?  Find the m's and find the l.  So you control the exit type of error, essentially that is related to the indicator function of  the path going out.  That will produce an O1 error.  So essentially that's why you have to compute to control this via the exit probability.  Then the weak error that is more or less what we have discussed before with the diffusions  and then the statistical error that doesn't change from the discussion that I made before.  So the issue of this is here is that you have to do things correctly to simulate couple  pairs that have the right law and are correlated in the right way.  And of course being able to accurately estimate the different quantities in the multilevel  construction.  This is particularly nonstandard here from respect to the diffusions that I discussed  before for the following reason.  The paths, right, they're taking values on a lattice.  So when things are going to converge, right, I mean as you discretize the L, the delta  t or increase the L if you wish, the paths are going to be closer and closer.  But you know, in diffusions closer and closer may mean you're in a continuous space, so  they may mean actually getting closer and closer.  But here you have a lattice.  So what does it mean to get close in a lattice?  You're far or you're there?  So actually closer in this setting, it means that the probability of going far is going  to zero, right, and the probability of being equal is getting larger and larger.  What's the problem with that?  That the problem is that this creates immediately high kurtosis.  And high kurtosis is bad news in the multilevel Monte Carlo context because it means the estimate  of variance is difficult.  So one thing that appears here naturally as opposed to the diffusion case that I discussed  before is you have to work a little bit farther to estimate the right variances in your system.  Remember the m's are driven by the variances.  So you need good estimates of these quantities.  Otherwise your estimates on the m's are wrong.  And if the estimates of the m's are wrong, bye-bye Monte Carlo, right?  So that's important from the practical point of perspective.  I think it's in delta and the index and the...  Little n, yes?  It's in delta that you use to balance the probability of s.  Yes, everything goes there.  So you know that.  Everything goes into there.  I mean, each level in principle here contains a delta Tl but contains also a delta L.  Because remember, there was a...  There was this probability of going out, right?  So you may prescribe this delta L as a function of the L and also a delta L maximum per level.  So this thing goes as before like delta T0 divided by 2 to the L, for example.  But this doesn't go that fast.  This is chosen in a different way.  So these two things combine to give the method.  Is that clear?  So these two are varying.  All right?  Thank you.  Good.  Thank you.  Okay.  So this is an example, the example that I showed at the beginning.  And here you see that from the exact method, essentially, which is SSA, the error bound  versus the tolerance essentially goes like total to minus 2.  So essentially, this is the complexity that we have.  And if you use this hybrid multilevel, you can still keep the same complexity and what  you gain is a constant but the constant can be substantial.  All right?  In this case, it's 10 times faster.  And this is actually the usual plot that you would see in our papers.  It's the control of the weak error.  So what the user gives in this context is a tolerance but he also prescribes the confidence  level.  So in this case, the confidence level was prescribed to be 5%.  So when you do your simulations, essentially, you expect to see around 5% of them going  above the prescribed tolerance.  So you see here, the tolerance is getting lower and the corresponding sample error for  the different runs of the material, Monte Carlo, distributed like this and then essentially  these numbers here are the percentage of those above the bound.  Okay?  So you change the quantity here, these exceedance values should change as well.  All right?  So you can actually have, depending on the problem, different type of gains.  I'm not going to discuss this too much but you can go above 10, of course.  You can get 10 to the 4, things like that.  But the plots look like this always.  You're improving the constant.  You're not improving the rate.  The rate was optimal to begin with.  Right?  So it is actually a pain to introduce a bias approximation with respect to the unbiased  method.  Okay?  And here's an example of our reaction splitting method where you actually do this differently  in the different reactions to get even farther.  Okay?  So this is, I think that's it.  So summarizing here, this is a different type of idea for multilevel adaptive, right?  You're switching between methods.  Okay?  You can also do the delta T adaptive as well but here you're switching between methods.  And you're going from an unbiased method into a bias method to improve the constant in the  complexity.  All right?  Any questions?  Okay.  So now we discuss discretization from another point of view.  We're going back to uniform discretization but now we ask a different question.  We say we have been driving our discretizations with this type of geometric type of refinements,  right?  I mean you get an initial discretization parameter, then you divide it by a given number, then  you divide it again by the same number, and so forth.  Okay?  This is natural within, again, numerical discretization, right?  You try with a mesh, then you start to divide it.  Try with another mesh.  It's not optimal.  It's not optimal.  That's a natural question to ask, right?  I mean, should we choose the HLs or the delta TLs like this or should we not?  Answer?  No.  It's not optimal.  Okay?  So with some assumptions you can show that this is not optimal.  This is the work we did in this optimization of hierarchies for Matias Monte Carlo.  And this other work that is here, this continuation of Matias Monte Carlo addresses this problem,  right?  That these variances that you need actually to drive and get the optimal ms are not known  a priori and somehow have to be found on the fly, right?  So you need to create an algorithm that not only is adapting the ms but is trying to find  these parameters that are needed to find ms, okay?  So that's what the continuation is doing.  And this is implemented in the lecture in this code, MIMC-LIVE, this library, MIMC-LIVE.  So let's go to the assumptions then to focus.  So what do we have here?  We have a Matias Monte Carlo method as before, but now the levels, GLs, are going to be driven  by these HL parameters of discretization and we do not say that they are geometrically  distributed a priori, okay?  So we're going to let them vary.  So we're going to take L capital plus one of those and we are going to assume that essentially  the weak rate is algebraically convergent.  The strong rate again is algebraically convergent.  And the work, for example, is algebraically divergent.  That's in the linear PDE that we have before.  And then what we have to find are not just the ms but also the Hs, okay?  Is that clear?  So as you know, the constants and the Qs here will depend on the problem and the solver,  as we discussed before.  And these are just a couple of examples, right?  In Euler-Maruyama you have Q1, Q2 equals one.  And if you have the other PDE example, we already discussed that the Q1 was two and the Q2 was  four.  Gamma again is defined by the solver and it's obvious in Euler-Maruyama, it's just one.  So I don't stop into that.  I just try to show you that these things can be found.  And if you have no theory at all, right, if you're desperate, you can just try to fit  something to the convergence as you see, right?  And or the explosion in the solver.  So you can actually measure things as well.  So two useful quantities here is the quotient between the strong convergence and the explosion  of the work that is chi and the eta is the quotient between the weak convergence rate  and the explosion of the work.  Of course, the best case for multivariate Monte Carlo is when chi is larger than one,  the limit case is chi is equal to one.  And the sort of algebraic improvement with respect to Monte Carlo is chi less than one.  All right?  And again, you see that if you change the solver, you may change from chi greater than  one into chi less than one as I described here.  In the PDE example, if you use an iterative solver, you get 1.34 and if you go into the  direct solver, you get 0.89 already.  So you can influence the convergence rate of multivariate just by changing the solver  like that.  Okay.  So the optimization, how does it look?  It looks like this work that was the work that you discussed before, subject to this constraint.  This is a bias constraint less than this fraction of the tolerance.  And this is the statistical constraint, which means that the total variance, right, is less  than theta times tau divided by c alpha.  And c alpha, again, is the confidence constant that is coming from the user, essentially,  when he imposes the confidence level alpha.  Okay.  So the good news is that you can solve this problem.  And of course, you can even verify the normality, how far you are or how close you are when  you do actual computations.  And here's a result somehow for chi equals one.  You actually have, is a case where you actually have geometric discretizations.  But I have to warn you, these geometric discretizations, you see, h a is somehow beta to the capital  minus l.  So it's a subdivision that we discussed before.  This doesn't depend on the levels, but it depends on the total.  So it's not your usual geometric discretization that you would think of.  Okay.  And the beta, actually, is not two or three.  I mean, it's something that may even have an influence from total.  Okay.  So the spacing between the levels is not trivial either.  All right.  This is one more thing that is interesting here.  The optimal choice of the splitting parameter, which is the fraction that you devote to the  bias versus the fraction that you devote to the statistical error, actually goes to zero,  sorry, goes to, oh, come on, goes to, I mean, the fraction of the bias goes to zero as a  total goes to zero.  Yes, that's right.  So essentially, if you see here, theta is the fraction that you give to the statistical  error and one minus theta to the bias.  So theta essentially goes to one, why, because theta is fixed and L goes to infinity.  So this theta is just going to one.  All right.  So asymptotically, you're putting most of your error, you know, your error allowance to the  statistical error, not to the bias.  And of course, the level behaves like log total and with this constant.  But for non-Ka equals one, things become less simple.  You can still find the relations on the theta and you can find bounds on the number of levels  with respect to the different parameters as this shows.  Okay.  Nothing really changes from the classical results in terms of exponents, but you can identify  how the input constants affect the actual complexity in the actual constants.  So you can play tricks to reduce these constants knowing that they exist and the shape they  have.  Okay.  So the shape they have in these constants in the different cases are like this.  So you can actually see that in some of them, for example, the variance of the level zero  plays a role, but in some others they don't.  So here in this setting, if you're here, you may get an advantage by doing a special trick  for level zero.  So whereas in the other cases, you may not.  Okay.  So like I said, the exponents are the same as you have in the classical theory, but now  you identify the constants as well.  And you identify how you have to space the levels too.  And that is more or less the discussion that I have here.  These are not going to be necessarily nested because you have this beta that behaves in  a strange way.  And the H0, also the first mesh, also behaves in an interesting fashion too.  Here's a discussion on the continuation of multi-level Monte Carlo.  And again, the plots always look like I discussed before.  You have a direct solver here and an interactive solver.  You have the prescribed error.  And here's it versus the tolerance.  And you're checking essentially that when you run over and over and over the multi-level  approximation, you get a fraction of those that go above in accordance to the more or  less the confidence level that you prescribed.  Now coming back to this question of finding the variances and other relevant parameters  to do optimization, in the continuation of multi-level Monte Carlo, what we do is to solve  a sequence of problems until you hit the problem that you want.  So essentially, what happens?  We know that in the best possible case, the explosion of the work in multi-level Monte  Carlo is told to be minus 2, right?  Best possible case.  It could be even worse.  And you can have a long talk here.  What that mean?  The worst is like this.  So if you actually solve now a sequence of problems where you first begin with told zero  you'll say, and then you go times told zero divided by a certain fraction, right?  Call it F, and then keep going, right?  Until you hit the actual tolerance you want, you can sum these works here, and you can  get a choice on this fraction so that the actual sum of these sort of preliminary works  does not create a huge overhead on the last one that you want to solve.  That's the first idea.  And then what happens here?  That when you start solving these sort of axiomary problems, you start learning the  constants that you need, okay?  So somehow by creating these fractions, you do not overshadow the actual work that you  want to do optimally, but you start sort of uncovering a veil over the variances that  you need and the other things that you may need to find the optimal elements.  That's essentially what is behind the continuation of Monte Carlo, sort of a motto with respect  to a tall type of idea.  Not really complicated, but it has to work quite nicely.  Of course, on top of that, there are tricks to actually improve a little bit the variance  estimates on the deeper levels, right?  What happens statistically is that in the course levels you have a lot of samples, just  by construction, right?  And Monte Carlo produces a lot of samples in course levels and very few in the deep ones.  But the deep ones are the really costly ones.  So if you make a mistake in the variance there, you make a mistake on the M and things go  really rotten.  So to get a better hold on the deep levels, you actually do some kind of Bayesian type  of estimate where you use as a prior model an extrapolation that is coming out of the  course levels and you enrich it with a few samples that you have from the deep levels.  That's what is behind too.  So it's not just one idea, but I mean, I'm just telling you the couple that I consider  most relevant.  So the good news again is that if you do all this sort of little learning process, the complexity  that you get is not much of an overhead with respect to the optimal line here.  That will mean that for this sort of pink magenta type of a curve, you are just writing the  cost that you would incur if you would know all these nice constants beforehand with respect  to finding them on the fly.  So it's not too much that you have to pay, really.  So that's with respect to the optimal splitting theta and it's not really much relevant.  So what is the second part of conclusions?  We didn't discuss much of normality, but we have a result using the linear failure in  this context.  Better splitting between Bayesian statistical error can improve the work a little bit and  this continuation of Monte Carlo tries to provide a more stable algorithm that allows  you to find the constant in between.  Geometric hierarchies are nearly optimal, but they are not the geometric hierarchies  that you are used to work with.  They come from this optimization.  They have specific values that may depend on toll as well.  So for ETOs, I included in the slide some assumptions and some results with respect  to the densities that I don't have time to comment on, but I'm going to share the slides  so if you have questions you can come to me.  Again I emphasize here that these files actually depend on derivatives of the final solution  with respect to initial conditions.  So this is what really makes the method non-adapted to the natural filtration.  So let's say this is multilevel in that context.  This is what I meant by strong convergence and the accuracy result.  This is what I meant by the work optimality result and I included a little bit more in  terms of complexity.  But OK, CLT2.  So last thing before I run out of time, a little bit of an activity with respect to  M and I just touch upon the classical expected value computation using just Monte Carlo just  to close it up.  So what's the point here?  You want to compute the expected value of X using IID sample averages and you want to  prescribe again this delta which is the confidence level and you want to prescribe also the toll.  So what is the point here?  The point is that usually you do this approximation based on the central limit theorem.  That rule, that will be a rule of the following type.  Essentially you are going to prescribe toll and you are going to prescribe delta and then  a number of samples M will come out.  It may not be so reliable when you have X that has low decay and density, it's essentially  heavy tails.  So for that reason we provided a rule that is a little bit more robust and essentially  combines the central limit approximation with other type of bounds.  Let me walk you through a little bit before I finish.  I have 10 minutes, right?  Yes, okay.  So this is going to be a little bit intense then.  So this is a general stopping rule, right?  Based on the samples you are going to estimate the probability of failure and based on that  we are going to accept or not.  So the probability of failure actually is the probability of your approximation not  being close enough to toll with the prescribed probability.  So this event of failure having a value larger than delta essentially.  So again as we described with the central limit theorem before the use of variance constraint  you can actually, if you believe that your sample average, your sample variance is actually  accurate enough you put a rule of the central limit type substituting sigma here which is  a true standard deviation by its approximation and if this actually is accepted then you  stop and you give your approximation.  These may have issues, okay?  So consider the following Pareto distribution.  The density is decaying like X to the minus 4.1, right?  For X larger than 1.  And you see here the behavior of this quotient, this quotient is the probability of the actual  failure versus delta the one that you wanted to impose.  So this quotient is larger than 1, this is bad news, this quotient is less than 1 or  equal than 1, it's okay.  So you see here tolerance versus delta, right?  So here it's like having more and more confidence in your computation and here is more and more  accuracy as you go to the left, right?  So as you let tolerance go to zero essentially for delta fixed things are going to be fine  but in the regime where delta is say a few percent and sorry, when total is few percent  and delta is sort of stringent and that's natural for engineering computations this  is not going to be nice, it's going to be actually in the red and it's not going to  give you the prescribed tolerance.  So there are many results in the literature where you fix the delta and you let total  go to zero.  What I claim is that essentially those results overlook this regime which is important in  the applications.  Okay, so the whole point here is to penalize this fact because we really don't believe sometimes  that we have such an accurate estimation of the variance to uncover this what we do essentially  is to use very same type of bounds that are essentially sort of better control rates on  the approximation in the central limit theorem and combining those with the sort of H-worth  type of expansion which is something that you would use when you have full confidence  that the tails are going fast essentially to zero then you get sort of an improved rule  and when you go above this rule you actually see in this example that you can achieve the  right prescribed failure rate with almost no change, almost no change you see in the computational  cost and the cost here is just determinants of samples because we are working with just  Monte Carlo in this context.  And this can be done also for another kind of not so nicely behaving random variable  and you see again that for small delta and relatively few percent total this is a bad  case but you can fix it using this algorithm and without incurring too much in an increase  of the work.  And finally this is an example where we actually expect the usual approach of central limit  theorem to work nicely and it does actually the work here is such that the failure probability  is below the delta but this algorithm also achieved that with that an increase, a sensible  increase in the work.  So in the nice cases this behaves well, in the worst cases or the bad sort of cases this  will have better control on what is going on.  Okay so that's a summary and I think I get to my end that's remarkable actually.  So thank you.  There any questions or comments?  Yes, Amar.  So I would be interested to know if we could apply this criterion for example in the case  where X is related to Rar event.  So that is a Bernoulli random variable typically which is 1 and 0 but 1 with very small probability.  What?  That's not the case we looked into.  I think that needs a measure change.  I don't think one can apply it directly here.  Just using some simple averaging of your output.  So our robust is the criterion to...  I mean what I'm saying is that it can be used but provided the first you do a measure change.  If you don't do a measure change all the estimates are going to be wrong, right?  You mean wrong inaccurate in fact?  Absolutely.  The whole thing is completely reliable, right?  So Rar event asks for a measure change to start with and then if the measure change is  sound, the important sampling is sound then you can use things like this.  But I wouldn't recommend using this in a Rar event without looking into an important sampling  algorithm first.  It can be combined with the same...  That's what I'm saying.  But I wouldn't recommend it to use it without the important sampling step first.  I have a second question perhaps.  Because again that regime is different, right?  I mean that regime is like for example 0, 1 with a very tiny probability of achieving  1.  I mean here it's like unbounded, right?  You can get unbounded variables with a certain probability and this is the case, not that  fast.  So you can achieve them with the probability that it's not so negligible and this is what  gives you trouble in the estimation of the basis.  It's not the same regime as the...  Have you looked in the direction of concentration of measure inequalities like we know...  Well, that's also not applicable here in the sense that for concentration of measures at  least if you use classical large deviations you need exponential moments.  These guys do not have any, right?  I mean we use them in the Cherno for the Taulib but here you don't have an exponential moment.  So if you have a concentration it's too little.  A short comment.  Maybe the applicability of your method is connected with a kind of ratio between the  expectation and the variance.  Because the case of Fourier events is that 1 over square root of p goes to infinity and  p goes to zero, something like that.  The relative value of the expectation...  Yes, yes, but I mean that's also what defines the regime of red events, right?  Yeah.  But in this case again I emphasize that is the tail that is creating the issue and we  create sort of something that works more or less uniformly with respect to that and does  not address at all the question of the red event.  The red event should be addressed in a complementary way.  At least it's my humble opinion.  Other questions, yes?  Have you looked or is there any chance to extend this framework to the case of PySYS  adaptivity?  I'm thinking of SPDE where you would like to have a discretization that has some local  adaptation in space.  But this is not about this last part, it's the previous part.  Yes.  Well, I'm kidding.  Okay, okay, okay.  This is about the...  Again, the adaptivity would respect to...  Pathwise.  Pathwise, okay.  So what we discussed was pathwise but in time.  Yeah, but the time step is fixed.  No.  Independent of the...  No, no, no, no, no.  The time step is delta T of T and omega.  Delta T is fully stochastic.  So it cannot have two features that come with omega, pathways.  So yeah, that's actually the most general but you only use it if you're obliged to.  That's my advice.  That if you can use, for example, the terministic type of activity, that is cheaper.  But you need some prior information about the problem that you're solving.  Right?  Okay, thank you.  Did you have any photo questions?  If not, so let's thank again the speaker.  Thank you.  Thank you.  Thank you.  Thank you.  Thank you.  Thank you.
2017-07-20_Tempone_part1_0_47.png 2017-07-20_Tempone_part1_48_95.png 2017-07-20_Tempone_part1_96_155.png 2017-07-20_Tempone_part1_156_203.png 2017-07-20_Tempone_part1_204_275.png 2017-07-20_Tempone_part1_276_407.png 2017-07-20_Tempone_part1_408_1355.png 2017-07-20_Tempone_part1_1356_1451.png 2017-07-20_Tempone_part1_1452_1799.png 2017-07-20_Tempone_part1_1800_2099.png 2017-07-20_Tempone_part1_2100_2723.png 2017-07-20_Tempone_part1_2724_3191.png 2017-07-20_Tempone_part1_3192_7235.png 2017-07-20_Tempone_part1_7236_8219.png 2017-07-20_Tempone_part1_8220_19403.png 2017-07-20_Tempone_part1_19404_20063.png 2017-07-20_Tempone_part1_20064_25439.png 2017-07-20_Tempone_part1_25440_36191.png 2017-07-20_Tempone_part1_36192_37319.png 2017-07-20_Tempone_part1_37320_37511.png 2017-07-20_Tempone_part1_37512_41327.png 2017-07-20_Tempone_part1_41328_41699.png 2017-07-20_Tempone_part1_41700_52751.png 2017-07-20_Tempone_part1_52752_53123.png 2017-07-20_Tempone_part1_53124_55787.png 2017-07-20_Tempone_part1_55788_61199.png 2017-07-20_Tempone_part1_61200_63695.png 2017-07-20_Tempone_part1_63696_63851.png 2017-07-20_Tempone_part1_63852_64343.png 2017-07-20_Tempone_part1_64344_69779.png 2017-07-20_Tempone_part1_69780_72659.png 2017-07-20_Tempone_part1_72660_73235.png 2017-07-20_Tempone_part1_73236_73355.png 2017-07-20_Tempone_part1_73356_73595.png 2017-07-20_Tempone_part1_73596_73787.png 2017-07-20_Tempone_part1_73788_74195.png 2017-07-20_Tempone_part1_74196_74399.png 2017-07-20_Tempone_part1_74400_83015.png 2017-07-20_Tempone_part1_83016_83243.png 2017-07-20_Tempone_part1_83244_83723.png 2017-07-20_Tempone_part1_83724_84215.png 2017-07-20_Tempone_part1_84216_84299.png 2017-07-20_Tempone_part1_84300_84455.png 2017-07-20_Tempone_part1_84456_84995.png 2017-07-20_Tempone_part1_84996_85871.png 2017-07-20_Tempone_part1_85872_86519.png 2017-07-20_Tempone_part1_86520_92579.png 2017-07-20_Tempone_part1_92580_93503.png 2017-07-20_Tempone_part1_93504_97667.png 2017-07-20_Tempone_part1_97668_99419.png 2017-07-20_Tempone_part1_99420_104123.png 2017-07-20_Tempone_part1_104124_104879.png 2017-07-20_Tempone_part1_104880_109763.png 2017-07-20_Tempone_part1_109764_113483.png 2017-07-20_Tempone_part1_113484_114803.png 2017-07-20_Tempone_part1_114804_114851.png 2017-07-20_Tempone_part1_114852_114983.png 2017-07-20_Tempone_part1_114984_115139.png 2017-07-20_Tempone_part1_115140_115187.png 2017-07-20_Tempone_part1_115188_115331.png 2017-07-20_Tempone_part1_115332_116027.png 2017-07-20_Tempone_part1_116028_116255.png 2017-07-20_Tempone_part1_116256_116543.png 2017-07-20_Tempone_part1_116544_116879.png 2017-07-20_Tempone_part1_116880_121379.png 2017-07-20_Tempone_part1_121380_121523.png 2017-07-20_Tempone_part1_121524_121919.png 2017-07-20_Tempone_part1_121920_122015.png 2017-07-20_Tempone_part1_122016_122087.png 2017-07-20_Tempone_part1_122088_122183.png 2017-07-20_Tempone_part1_122184_122351.png 2017-07-20_Tempone_part1_122352_122939.png 2017-07-20_Tempone_part1_122940_123119.png 2017-07-20_Tempone_part1_123120_123179.png 2017-07-20_Tempone_part1_123180_123491.png 2017-07-20_Tempone_part1_123492_127307.png 2017-07-20_Tempone_part1_127308_134171.png 2017-07-20_Tempone_part1_134172_134411.png 2017-07-20_Tempone_part1_134412_134483.png 2017-07-20_Tempone_part1_134484_134615.png 2017-07-20_Tempone_part1_134616_136259.png 2017-07-20_Tempone_part1_136260_136595.png 2017-07-20_Tempone_part1_136596_137723.png 2017-07-20_Tempone_part1_137724_138419.png 2017-07-20_Tempone_part1_138420_140255.png 2017-07-20_Tempone_part1_140256_140615.png 2017-07-20_Tempone_part1_140616_141731.png 2017-07-20_Tempone_part1_141732_142187.png 2017-07-20_Tempone_part1_142188_143039.png 2017-07-20_Tempone_part1_143040_146963.png 2017-07-20_Tempone_part1_146964_148871.png 2017-07-20_Tempone_part1_148872_151031.png 2017-07-20_Tempone_part1_151032_153839.png 2017-07-20_Tempone_part1_153840_154847.png 2017-07-20_Tempone_part1_154848_156743.png 2017-07-20_Tempone_part1_156744_157067.png 2017-07-20_Tempone_part1_157068_158843.png 2017-07-20_Tempone_part1_158844_159059.png 2017-07-20_Tempone_part1_159060_159107.png 2017-07-20_Tempone_part1_159108_159155.png 2017-07-20_Tempone_part1_159156_159225.png
Thank you very much for coming with me today. And I'll give two lectures. The first one will be related mainly to multilayer Monte Carlo and the second one we will touch upon multi in the expressions of the things that I'm going to say today, but not in an adaptive fashion as we will see. So first go then to discuss a little bit about the activity within multi-level Monte Carlo. Let's see how this goes. So the plan is as follows. I will, well, these are some pictures from the place I come. This is a chaos actually. It's a King Abdullah University of Science and Technology. Some of the students are here today. And let's see. This is a plan. So most of you are very well acquainted with Monte Carlo and possibly multilayer Monte Carlo as well. So I will recap a little bit just to set up the notation at the beginning. So we will speak the same language. And then we go into different versions of our activity. As you would see, there are many ideas for adaptive approximations in multilayer Monte Carlo that we can see. And for that reason also, I will touch upon different models. I mean, one of them will be computations with random PDEs, that is PDEs that have random coefficients. And the second one will be diffusions. These are etostochastic differential equations. And the third one will be pure jam processes, say stochastic reaction networks, whose trajectories are piecewise constant. So I think this will motivate a little bit different ideas that we will have in this. And finally, we will discuss a little bit of the contribution multilevel and optimal hierarchies. Hopefully, if the time allows, I will also discuss optimal stopping with non-asymptotic rules within Monte Carlo. So I'll try to really be adaptive in this process. So let's begin with our first example. This is a linear elliptic PDE that is very popular in the context of random PDEs as a model case. So here, the coefficients a and f are random processes, stochastic random fields. And mainly, we'll concentrate on the effect of a, because a is with the two, right, the one that acts non-linearly on the solution u. And out of the solution u, we are interested in computing a statistic, which is real value. Again, as an example, in this lecture, I'm going to go through the computation of the expected value of a functional of the solution here depicted by psi of u. Psi could be linear or non-linear. We have some restrictions on smoothness on psi. But for the moment, think of a linear functional that is the magnetic acting on u. Again, this is not the only thing that you can compute with Monte Carlo, but this isn't the one that we will be focusing mostly through the lectures. All right. So the goal is set. We want to compute this expectation. Now this expectation involves the solution of a differential equation. Seldom, do we have the actual exact solution of u of our disposal? So instead of u, we're going to approximate it by uh, say, where uh comes through a discretization of this equation. In this case, of course, it will be only space discretization. But in general, you may have time as well or even other parameters. So think that you have prescribed the use of your favorite method. In this context, you could define it elements, find a difference, find a volume, so whatever. And that discretization is indexed by the parameter h. All right. So h will govern in two aspects of the solution. One will be the approximation characteristic of uh to u, but at the same time, the cost of evaluating uh. With the assumption, of course, of consistency, meaning that when h goes to zero, uh converges in the proper space to u. And at the same time, the cost of evaluating uh goes to infinity. Right? Otherwise, we will let h go to zero and get the exact solution. Now, since we were interested in this function actually, which is a real value random variable, we will approximate, I mean, the notation denote the approximation by gh, where now h is the discretization parameter. All right? So that's hopefully a natural notation. Let's go through the error then. The difference between what we would like to compute, which is the expected value of g minus what we actually compute, which is a sample average based on IID samples out of the inexact quantity, which is gh and not g. Right? So we split upon two main sources. Right? One, again, the trick is to somehow subtract the expected value of this approximate quantity to split the difference between this and this into two things. Right? One that we call the bias error, the other is a statistical error. And for the bias error, we have essentially an estimate of this kind, which is coming directly essentially from the deterministic analysis, right, point-wise, and then integrated upon this expectation. So essentially what we have is that if your method has order w here with respect to h, which is the discretization parameter, when essentially the constant that comes from the deterministic analysis will be random, of course, because it will be depending on a and an f. But once you take the expected value, you get the average of such a constant, which is no longer stochastic, it's just a deterministic constant. All right? So this analysis holds as long as there are no changes of order with respect to the randomness in the problem, which is a reasonable assumption. But if that happens, actually, this will not be the correct way to draw things. But again, this is an actually way to explain how this works. So let's again consider that this w doesn't change with respect to omega, and therefore we just get a constant here times h to the w. Okay, so this is more or less a standard discretization error. Second term is the standard statistical error, right? And this error essentially is mainly governed by the number of samples in. And I say mainly because the random variable that you're actually sampling depends on h, right? So this is an unbiased error, correct, because we take expected value here and then some sample of the same random variable in an ID fashion. So we only need to control its variance, right? So and that is motivated by the central limit here, essentially properly scaled, right? This will be approximated by a normal distributed random variable. And this allows to write then later the constraint of having this error smaller than a certain quantity in probability, you know, as a constraint in the variance, right? So as we can see, again, we have two contributions here, one coming from h, the other coming from m. And as long as the problem is a rational well-posed problem and there's nothing strange, these variants is uniformly bounded with respect to h, right? And therefore they claim that m is essentially the quantity that governs the statistical error, okay? And if the problem has a singularity of some kind, maybe this is not completely true, but I mean, again, in the examples that we're going to visit, this variance is going to be bounded uniformly with respect to h, okay? So two sources, one is governed entirely by h and converges algebraically with respect to h. The w here is determined by the numerical method that you use and the regularity of the problem that you're solving. And this m is very robustly coming into this analysis, right? We only need to ensure that this variance is bounded, essentially. So now that we have the two main errors identified, we try to see what is the corresponding complexity, which is the computational work that is needed to achieve a certain accuracy with Monte Carlo in this context for a prescribed confidence level, okay? And to do that, essentially what you do is to minimize the work subject to the total error constraint, okay? So the total error constraint, again, improbability, because this was essentially an analyzing probability, is the sum of the bias plus the statistical error. The bias was governed by this, the statistical error governed by that. And the total work, naturally, is what? Is the product between the number of samples that you have to solve for, right? Sometimes the expected cost per sample, okay? And again, this is an example, right? And the example that we have illustrated is a linear elliptic differential equation, right? And in that case, what we have to do is to solve a linear system every time we need to create a new sample of the random variable. So if v is the dimensionality of the space that you're solving, right? I mean, could be 3D, for example, 2D, 1D. Then h to the minus d, if h is a discretization parameter, is essentially the number of degrees of freedom that you have in your discretization. So h to the minus d is a degrees of freedom. And then to the power gamma is essentially the cost that one incurs by solving the system, okay? So gamma equals one will be an optimal solver, and gamma equals three will be the Gaussian elimination, okay? So gamma can raise between one and three. And of course, h to the minus d is the number of degrees of freedom in your system. Is that clear? Okay? So this adapts to, again, to the type of solver. If you use one type of solver or another, this gamma may change, okay? Let's see what else is to be said. Yes, okay. So I mentioned the following. I said that this cost was understood in the mean sense. And why do I say that? Well, if actually you have an iterative solver, right? And the spectrum of the operator that you're solving with is random, right? Because the a is the feasibility coefficient, was random itself, the number of iterations may depend on the realization, right? So what you're really interested in is the mean cost, not just the cost. I mean, the mean cost will be representative. And there will be variations with respect to each of the realizations. Is that clear? Okay? Now, if you use a direct solver, this will not be the case. But okay, I just want to emphasize a little bit on that. Now, if you actually do the minimization of this bound and subject to this constraint, you get this classical cost with respect to the tolerance that you're imposing. Okay? In the context of the central limit theorem, the confidence level is only acting on the C2. So it really plays a role in this constant. I mean, in this exponent, it's just acting on the constant. And that's why you don't see it. But again, everything now depends on the accuracy you want to achieve and can be split into two contributions. One that is told to the minus 2 and the other is told to the minus d gamma divided by w. Now the first part, told to the minus 2, is the classical cost in terms of number of samples of a Monte Carlo method. This is essentially what is needed to make this less and tall. So that will be the cost that one will have to pay if each sample cost me over 1. The second part of the cost is coming from the fact that I need to make h sufficiently small so the bias constraint is satisfied, meaning that this term is actually less and tall. Okay? So essentially what you see here is the number of samples times the cost per sample, which is essentially what you were seeing here as well. Okay? Now, what Monte Carlo tries to do is to go beyond this bound. So what we will improve somehow is this part and this will still be a barrier that we will not be able to cross at least with Monte Carlo methods. So just to give the notation again for multilevel Monte Carlo, let's recall this construction of Henry and Giles. Let's take a constant beta for subdivision. For all purposes, think of beta equals 2 for the moment. Okay? So what we say is the following. Let's just take in a single h that was meant to satisfy the bias constraint. What we're going to do now is to take a number of such h's. In principle, if you come from analysis, from numerical analysis, this is a natural way to proceed because every time you are actually trying to solve a differential equation, you try with a certain mesh. If you're not happy, you subdivide it right and you test the stopping criterion. And if you're not happy, you subdivide again until you meet the stopping criterion. Normally, what you do is that once you get to that level that you accept the error to be smaller than your prescribed accuracy, you throw all the other things, right? You don't care, essentially, and you just work with the last accepted guy. Here on the contrary, you're going to keep those for a certain purpose. And the certain purpose is to produce an effective control variant. So what one is doing here, again, let's go through the definitions first, is to introduce these delta operators that are differences between consecutive approximations. So essentially, this, remember, is a functional evaluated using a mesh of size hl. And this is a functional using a mesh of size hl minus 1. This is the difference between two random variables. And we are assuming, of course, that they are not independent. They are going to be sampled according to the same omega. That's fundamental in this construction, because that is what will make this difference actually go to zero, essentially, when l, bit l here goes to infinity. That again comes from consistency. In l2, this is what we're going to mainly need for this to be working fine. So let's go again to the construction. The construction of these hls is done via a geometric discretization, meaning that you start with a given mesh. If you're not happy, so divided by 2. If you're not happy, so divided by 2 again, by 2 and so forth. That creates this hl based on h0 divided by beta to the minus l with beta equals 2. Is that clear? So now that these deltas have been defined, then one rewrites. Remember this was the guy that we were trying to use as an approximation to this one, so that the bias between these two was less than tall. Now this is somehow the random variable that one is going to apply control variance to, and before doing so, one rewrites this expectation as this other sum. So this is an exact representation, and the sum is running overall levels of approximation from zero to the level l capital that satisfies the bias constraint. Yes? Now instead of using multicadload here, one instead uses multicadload independently for each of the terms coming into this sum. All right? So, and that's sort of a scary proposition because the first time you see here, this is saying at the beginning I only had one expectation to compute. Now I have l capital plus one. How is it going to work? Well, okay, before we look into that, let's see how this looks. And as you see, this, which is actually the estimator for multicadload, looks precisely like the sum of l capital plus one multicadload estimators, where each of them now is sampling something different, each of them is sampling these deltas. And of course, remember that for l equals zero, actually, this is not the difference, it's just g evaluated on the course's approximation. Okay, so you can think of this as, you know, g is h0, right, if you think of hierarchy approximation, and then the other things are just corrections, right, that are being added to the course's level. Now that we have this, how do we motivate, at least first intuitively, why is this going to work better than the original multicadload approach? Well, let's see, the bias is not going to be changed. We said, essentially, if this L is chosen so that, you know, this bias is less than told, what we are rewriting, essentially, is this quantity. So if you're going to actually improve this sampling is because we're going to pay less in the samples, essentially. The difference will not come from the bias. Yes, please. So how do we choose for each level? Oh, this is coming. This is coming. This is coming, don't worry. This is coming. Okay. Now, the variance here is the variance of multicadload, right? So what did we say that the constraint, you know, the statistical constraint looked like? Well, essentially, thanks to the construction, based on independent samples, the variance of this estimator is just one over ML capital of the variance of this GL capital. Now, you think that the GL capital is more or less close to G, so more or less you can expect that this looks like the variance of G divided by ML, right? So this is essentially saying that in the deepest level of approximation, the capital, the one that satisfies the bias constraint, you need to use total to minus two samples, right? And that is what you expect from multicadload. If you, on the contrary, use multi-level multicadload, right? The variance for the estimator, again, by the construction of independence, not only between the samples that we use on each level, but among levels, right? We get some of variances across levels. The first variance, the one that is coming from the level zero, actually, looks very much like this term, right? This is not a difference. This is not a difference as well. But the other terms are actually one over M, right, times variance of differences. And what's the point here? Remember, as little L gets larger and larger, we hope to have convergence, okay, just by consistency. Now, in which sense we expect to have convergence in this setting, what is most useful is to consider not only you have convergence in the mean, but also in L2, which means that these variances of differences are going to zero, right? This is where multi-level multicadload becomes effective. And if they actually become very, very small, it means that for really deep levels when L is large, right, each sample cost me a lot, but the variance is correspondingly small. So the number M that we have to use is really small. And this is where you get real advantages, because you are sampling much less in the most costly levels, all right? Now, you have to check some conditions usually to make sure that the choice of the level zero say sound, but again, the whole thing can be put in terms of these inequalities that again indicate that L2 is very relevant in the choice of the multi-level parameters. So far so good? Any questions? You said there was to be correlation between the different levels. No, there's not. Actually, no, they are constructed totally independent. Okay? That's important. Yes. From the independence in the construction, you get immediately this sound. Okay? Otherwise, you will have extra terms, all right? It's not that it's impossible to do that, but the construction here is done on independence in some chance. Okay. So cartoon type of approximation essentially, we are somehow moving samples from sort of deep refined levels that are costly into levels that are cheaper and less refined. Okay? Okay. So if you want to go beyond that intuition, you need some assumptions, right? And these assumptions may or may not hold in your problem. You have to verify that they are satisfied and sometimes you may have precisely this kind of convergence and some other type, I mean, depending on what problems you have. So here we are again going to concentrate on this algebraic type of convergence with respect to H. Remember, H was the parameter of discretization. You can think of delta x, delta z, delta y in your linear elliptic PDE. And we were assuming that the HL was essentially H0 times beta to the minus L, right? A geometric subdivision. So this first equation here is nothing different from what I wrote before, right? It means that the error in the mean is H to the w, essentially. You see? So nothing is different. It's just a rewriting of what we have before. Now again, this w, which is the exponent of the weak convergence, depends on the assumptions on your problem and the numerical discretization that you're using and has to be found for each problem, right? The variance assumption again is indicating some kind of error rate, but this error now is measured in L2. And instead of taking a difference between g and gl, now we are taking, remember, the difference between gl and gl minus 1 when we put it inside of the difference. And we have a strong rate of convergence that is s, which is not going to be w, right? In general. So for the last one, we have the cost, the mean cost of evaluating a difference on the level L, right? So this is delta gl. And this again entails the cost of evaluating gl and gl minus 1 because we're going to take the difference. Usually when assumptions that gl is the one that dominates, so this goes like that. And this is nothing from what we said before again, because before we said that the cost was H to the minus d gamma. And again, this is the same thing. Just rewritten in terms of this beta and this L. I emphasize again that the constant c, w, and s also depend on the problem and have to be found. And remember, gamma may depend also on the solver that you use, not just on the problem and the discretization. So if you change the solver, gamma will change even though you keep the rest of the things constant. As an example, for a smooth linear PDE of the type, linear type that we discussed today, if you use multi-linear piecewise continuous finite elements, we will have 2w equals s equals 4 and of course, depending on the solver, use gamma ranging from 1 to 3. So what's the total work? And here we start getting close to your question about the m's. The total work will be found in the same way as we did with Monte Carlo. We have L capital plus 1 levels. On each level, we have to use ML samples and the cost per sample. It's the work of delta GL in a expected sense. So what do we do? Well, remember, we have to find ML. So if you fix L, essentially, you can write an optimization problem and I'm going to go into that a little bit more when I go into multi-index Monte Carlo, but we can just say it again. I mean, remember, if you fix L, the bias constraint is not going to be important and what you want to minimize is minimize with respect to the ML's from L equals 0 to L capital. What? The total work. And this is the sum from L equals 0. I hope it is possible to see this. There's no reflection or anything like that. So L capital ML and this work delta GL. And then, of course, there's a constraint. The constraint is the statistical constraint. The statistical constraint, again, is the variance being less than the total square. So what is the variance again is just 1 over M0 say variance delta, actually G0. Just write it differently. Plus the sum from L equals 1 to capital of the variance delta GL. And this is less than total square. Of course, there's a constant here in between, but I'm not writing it. The constant comes from the confidence level that you want to impose. Okay. So the control variables here are ML, capital plus 1, right? And here you see how they act. So you can just write the last of relaxation and you realize that you can find these ML's. And ML's actually are going to be, at least if you make a real value relaxation first, this will be total to the minus 2. And then there's a constant here that is coming out of the levels. And this essentially is the variance of delta GL, I think. You may correct me if I'm doing this wrong. And then this is the cost, essentially. And this is for root here, if I'm not mistaken. Okay. So you can first write this lagrange of relaxation, solve it, and then say round it above, I mean, take the ceiling of this number to make sure that everything goes well. Okay. And from there, you get to that formula, essentially, assuming that the cost is not trivial, essentially, that you're not getting less than one realization per level. But that is a technical assumption. You can stay with this for the moment. So what do you see from here? You see that after you substitute the optimal number of samples, the cost only depends on L capital. Because these are problem, say, and discretization-based constants. So once you give the discretization method and the problem, of course, this work L and the variance L, essentially, are given. Maybe you don't know them, but precisely a priority, but they're given. And this is a sum of positive terms. So first thing that you observe is that when the tolerance goes to zero, L has to go to infinity. Right? Is it clear why? You have to satisfy the bias constraint, right? So the bias constraint, essentially, was h to the w less than toll. And L is coupled directly with this h to the w. So L, essentially, goes to infinity like a logarithm of toll, absolute value, say, that's a ballpark. You may have log log terms on top of that, but really it's log toll. So this is going to infinity as toll goes to zero. So this sum can only do two things, either stay bounded as toll goes to zero or diverge. If it is diverging, it's because this work is somehow growing faster than what the variance is being able to converge, right? But if it is actually finite, which is essentially the case where this variance is dominating the work somehow, we have a very interesting case. But it's sort of not intuitive the first time you see it as well. It's the fact that you get a cost that is essentially all of toll to the minus two. And if you remember Monte Carlo, when we discussed it first, toll to the minus two was the optimal cost with Monte Carlo when each of the samples, essentially, in the mean cost me all one. So essentially, this means that you don't feel the cost of discretization even though you're refining the levels, right, infinitely much as toll goes to zero. Is that okay? In principle, not, right? So somehow it means that even though we still have to sample of those at the level L capital, we only sample so few, so little samples so that that cost really does not influence the rest of the cost that come from the other levels. Because we still have to sample in that way. That is clear, right? Okay, so now that you see what is coming, let's see the result with these assumptions in mind. L capital, as I told you, is determined by the bias constraint and essentially is linked to toll by a log factor, right, where the log essentially comes connected with a weak rate of convergence and the logarithm of the subdivision multiplicative constant beta. And the samples ML are chosen optimally like we discussed before just by solving that Lagrangian relaxation problem. And after substitution, essentially we get this fantastic case, toll to the minus two, where S is larger than the gamma, right? So essentially both the work and the variance are behaving with respect to L in an exponential fashion, right? So the work is growing with the D gamma rate and the variance is decreasing with the S rate. So as long as S is larger than the gamma, this sum is convergent. And this is what explains the first result. The second is a limit case where essentially they outmaneuver each other. So essentially they balance each other. You get something that is constant here. So the total result here is L on this sum, L times a constant, and then you get L squared and L was log toll. So that explains the second line here. And in the last one, well, the work is growing faster than the variance is converging and instead of getting something that is convergent, you get some algebraic deterioration. But still this algebraic deterioration is better than what you would get with Monte Carlo. So you're decreasing from D gamma here in the explosion to D gamma minus S. OK? So you're still exploiting a little bit of the strong convergence. I mean, the convergence in it too that one assumed. So OK? So far for the introduction with the uniform methods. Now the question that we want to more or less touch upon is how to extend this Monte Carlo into non-uniform adaptive discretization settings and whether it's worth it. OK? And there are several ideas to introduce non-uniform discretizations or stochastic levels, if you wish. So let's begin with the... OK, well, this is... I put just to include in this slides and I already gave up trying because, I mean, the pages are not enough actually to describe so many things. And I rather just direct people to look into mics. Himongu said to keep up with the growth of the community. And this is actually what you get the best references in the... List of references in the area, so I'm not going to talk too much into that. So let's go to the activity then. And in this context, again, we're trying to approximate the expected value of an output quantity. In this case, I'm switching the example. I'm going to stochastic differential equations, itto differential equations. So in this setting, each of the paths evolved in finite dimensions with respect to time, right? So you can do this with the sum of the paths to a final value and then in that final value, you may evaluate the function g here that comes in this example, like a put option in finance and then you have to integrate this g versus the PDF of the final time. OK, so you can do this with Monte Carlo just by sampling many, many paths and again averaging as we discussed before. So here are some references again. I'm not going to stop and I'm going directly into the problem formulation. Here is an SE. This SE, if you have not seen it before, it's not a problem because you can think of this discretization via Euler-Marujama that is coming just in a second. So x is something that is finite dimensional now, it's not a function as in the first example. It evolves in RD. And g is a given function so that our goal is to compute the expected value of g of xt and so the evaluation at the final time of the solution of an SE. Again we are given an accuracy constraint tall and a prescribed confidence level and everything is driven by k-dimensional venerable process. All right, this has plenty of applications that I will not stop upon and what we're going to do is through this adaptive Monte Carlo and actually see that there are several versions now of adaptivity. The first one we did was with respect to the weak error but later we did something that controls the strong error and yeah, let's go into that later. So first discretization algorithm. This is just Euler-Marujama so if you have not seen SEs before, this is a new term that comes in. So remember it's coming through these. If you just think of the first part, it's just forward Euler right based on this function, the drift function a and the second is just the same type of discretization. You see you have to evaluate b on the left point, essentially the current point x bar but you multiply times these random variables. These random variables are normal distributed with mean zero and the variance is just the delta t that you're taking. So the delta t does not need to be uniform at all but yeah and again weekly you can just scale them to normal zero one random variables and the multicolline estimator looks exactly as we discussed before. It's an average of IID samples but the IID samples are not from the exact process but from the solution of coming out from this forward Euler. Forward Euler of course is not the only discretization that one can do in this context but it's the one that I'm going to exemplify. Okay so as we discussed before this can be split into tolerance for time discretization and tolerance for statistical error and this is a particular case as we discussed before the cost per sample here is essentially proportional to the number of time steps. The number of time steps essentially is proportional to one over a ton. This coming from the weak convergence rate one of Euler which is the same as Euler-Maluchin in this setting and this brings us to the complexity told to the minus three for Monte Carlo. Okay okay and now you want to introduce adaptivity. So what is adaptivity in this context and this context adaptivity means the ability to choose the time steps non-uniformly and this non-uniformly can depend also on omega. So the time steps such as t0 is 0 and the final time step is t capital are not going to be equally spaces and they may adapt to omega because the dynamics of the solution may depend on omega in a dramatic way maybe. So in this context then we will create realizations that depend on x bar, I mean here denoted by x bar and depend on this mentions delta t of omega. So again why do we do this? Well sometimes there are non-smoothness on A or B that ask to, I mean to deteriorate the usual behavior of the uniform method so ask essentially for the adaptivity in a natural way and what do we need to be able to produce such delta t? Well you need, at least with our approach, you need to rewrite the weak error in this other way essentially, you need to introduce an error density which is a function rho w that does essentially not depend on delta t so that the error, the weak error can be split essentially in contributions that are coming out of each of the different discretization steps. You see here if you discretize over the time steps, this is a sum over time steps of essentially delta t n square times rho w. Now the fact that delta t here is inside of the expectation is allowing for the choice of stochastic time steps because rho w may still vary dramatically with respect to omega, that's the point. If this rho w will not be varying dramatically with respect to omega then it's cheaper actually to take a non-uniform but deterministic discretization with respect to delta t, is that clear? There's another head to do this fully stochastic so whenever possible you avoid that. Now the funny thing is that in this context there's not only a weak density but there's also a strong density. You can actually make an approximation and you see here the power is one but here the power is two, you can actually write the L2 error here in the same fashion for Euler where now this rho s is not exactly this rho w, this is still something that you can approximate with computable quantities and it's able to distribute the strong error across the different intervals that come into discretization. In principle this is not obvious because when you square this you may get a lot of cross terms. So to be able to tell that actually the diagonal terms are the one dominant this needs some understanding. Okay, but what is the whole point now is that once you actually produce these error densities you can produce single-level adaptive algorithms. So single-level adaptive algorithms will try to find these time steps so that say for example a bias error is satisfied or even say this type of L2 error is satisfied. And what does it allow then? It allows for the following construction. If you want to create a multi-level adaptive method what you derive essentially the discretization with is with a sequence of tolerances. So if you have a reliable single level adaptive method the only step that you need to go from there to a multi-level adaptive method is just to create a sequence of tolerances and to say that the first level will be driven by what will be defined by discretization of my single level adaptive method so that they reach that accuracy. The next level is again refinements on those that achieve the next accuracy constraint and so forth for all the levels. So you see now each of the levels does not have a prescribed H. This notion is gone. Each of the realizations may use a completely different mesh that adapts to the characteristics of the omega that you have in mind. But still one will be able to extract the good properties of the variance convergence, the mean convergence and the rates of the explosion of the work with respect to the level. Okay? So what is optimal at least to show results on the algorithm in this context is to refine according and to stop according to these simple rules. Essentially what you do is to compute solutions and error indicators are n for each time step. And the error indicators are n in the context of this error density defined essentially as rho n times delta tn squared. I mean and rho n could be the weak error if you're controlling the weak error of the strong, if you're using the strong approximation that's not important. But again now this is defining the contribution from times the n in a given realization to the total work. And what we do to stop the discretization is to verify whether the maximum of these guys is or not above this. And this you see is the tolerance that we prescribe for the bias constraint which is what we want to achieve divided by the expected value of a number of time steps. So this essentially is the mean error that we would like to achieve on each interval. That turns to be optimal essentially to equidistribute the error along intervals. One caveat though, this n is not known a priori, right? Because this is the number of steps that you need to solve. So these have to be found on the fly. But this is a constant and this is another constant that are essentially probably independent. So what you have to remember here is that this is going to be not the same as this C of refinement. The refinement criterion is a little bit more stringent than the acceptance. So essentially if you see that you're not going to stop because some of them are going to be larger than that, well then you are going to refine many more than just those that not satisfy this constraint. Now another technical point in this construction is the following, if you have a refinement and has not been accepted and you mark some intervals in time for refinement, you need to produce essentially the corresponding interpolation values for the delta W's that drive your discretization. And these, you know, not to introduce any biases, have to be consistent with the points that you already sampled. And this is done according to Branian Bridges, right? That is the correct way to do this discretization. So this plot essentially shows that, I mean, that for a given level, for example, you use the blue discretization, right? And this entails evaluating W at all those points. But now you need to evaluate W in points in between. And of course, they are not going to be linear interpolants of the previous values, but they will differ upon. And this difference is actually given by the Branian Bridges, okay? So it's just conditional Gaussian sampling. That's nothing more. Okay, any questions? No? Okay, so what was the main idea again? Let the tolerance in the delta W algorithm define the hierarchy. That's essentially what you have to remember. And this hopes, this works in any problem. If you have a sound error density of PEs or EZEs, you can just use this approach and it will work. Like I said before, adaptivity is useful when you have issues with A or B or the process itself may have issues with respect to, say, boundaries where it is being, for example, stopped or reflected that those actually ask for adaptivity. You may have stiffness as well, and this is also nicely treated by adaptive algorithms. So again, if you have a stiff drift, you may choose at least delta T. If you don't want to go through this posterior global error control, if you just want to drive things by stability, you just may take delta T to make sure that you remain within the stability region of the method. And that again, if the stiffness is just coming from A, it can be done quite easily. So what is interesting here? The point is if only when you do a posterior error estimate, you show things about the discretization being efficient or not in the sense that the estimate for the error essentially gets close to the error you want to impose. But in this context of error densities, you can show actually, you show much more, you can show actually that the algorithm stops, that asymptotically you have normality, asymptotically you have accuracy. And again, what is asymptotic in this case is when tolerance goes to zero. There is no notion of H equal to zero because H is gone. So the only parameter that drives the discretization here is the accuracy constraint, and the accuracy is given by toll. So when you look at asymptotics in this context, it's when toll goes to zero. And of course, we have kept the confidence level constant in this discussion. We are not letting the confidence goes to one. So asymptotic normality, you also have asymptotic accuracy. I will just tell you later what this means. But it really means that when tolerance goes to zero, essentially the error is bounded by toll. And the complexity essentially is up to essentially a constant, you get the optimal complexity. So for the ethos to classic differential equations, this one discussion that I have not made, and I just put the warning here, it's actually working fine in the theory, but it's not obvious. The discretizations that you get are not adapted, I mean, they are adaptive, but they are not adapted. Adapted means that they are not adapted to the natural filtration of the problem, which is given by the linear process. Why? It's because the density, and I have not told you maybe, no, these densities essentially involve variations of the final value with respect to the current value. So essentially, they may need, for example, derivatives of respect to the position Xn of g of xt capital bar. This quantity is not adapted. It's looking to the future from, say, the current point n to final time. So in principle, you can discuss whether these things are going to converge to the correct solution. They do. But I mean, this is not true for any adaptive method that you can come up with. And no one has to be worried, I mean, at least careful about it. Okay, so that's why. So if you do examples here, for example, you have a drift singularity, this is just going to be simple in the case, in the sense that the stochasticity is really not relevant. I mean, here we have, say, deterministic time for the singularity, and the mesh is going to refine automatically around the value of alpha. We are not telling the algorithm that alpha is there, but it does find it. You see the mesh here is refined from the level, essentially, 10 to the 0 to 10 to the minus 10, just to control the error. Okay. But again, like I told you before, this is sort of the deterministic type of singularity. So in this case, it doesn't pay too much to let delta w to depend on omega. Delta w is not uniform, but it really doesn't depend on omega. So you can split the construction here in a first phase, where you actually find all these meshes that you're going to use in all the levels. And then you just run in the second phase, material multicolor increasing the number of samples that you need on each level until you hit the variance constraint. That is effective. Okay. And actually, what you can show is that the cost of the first phase to find actually the optimal discretization in Monte Carlo is much cheaper than the actual cost to solve the problem. So that's again, actually, when you have deterministic type of singularities that you can find the mesh very quickly. Okay. So in this case, you see an advantage from uniform to adaptive. You don't see much more of a difference between uniform Monte Carlo and adaptive Monte Carlo, because essentially what happens is that this singularity is not strong enough to modify the strong rate. It does modify the weak rate. But if you remember what I said before, the complexity in Monte Carlo is mainly governed by the strong rate. I mean, there's this independence of the weak, of course, but in the inequalities that I showed you before with these conversions of the sum, the strong rate was the one that was opposing the weak, the work growth, not the weak. Okay. So, but if you actually go to stronger singularities, so now we let P actually be beyond one half and, of course, you have to make sense that this is actually has a unique solution, blah, blah, blah, but it does actually, as long as P, I think is less than one strictly. So you can see here for P equals 0.5, that was the previous example, you don't see much of an advantage between multi-level Monte Carlo uniform and adaptive, but as you start raising the P, you clearly see the advantage between adaptivity and non-adaptivity. And again, I have to tell you one more thing, this has been done using the strong density, not the weak density. So this is the latest work that we have for adaptivity. It's actually on the strong, I mean, you use the strong density to rewrite this L2 error and you now create the levels according to this strong density. That's how this plot is being created. One more thing, okay, the alpha here was uniformly distributed between one fourth and three fourths, so actually the measures have to be stochastic, because the position of the singularities is being just jumping around. Okay, so I told you essentially how to do adaptive multi-level Monte Carlo approximation, yes please. Sorry, alpha is independent of w, that's right. So essentially you set alpha, you run the things, alpha is independent of w, correct. Yes, I mean, this was just a simple example. So people can actually reproduce and see what is going on. But yes, I mean, things can be more complicated because the singularity, as you're saying, can be coupled with the dynamics, okay. So it could be that actually the singularities appear because of some drift problem, right, with the stiffness, or it could be that x is getting close to a boundary and the boundary has an issue. So you can think of the boundary being penalized and penalization is drift stiffness, right, so yes. But that will not be able to be written in one line essentially, so yes, at least this motivates what you can gain. Thank you. Okay, any more questions at this point? I'd like to then summarize what we did on this topic. So I hope the idea of adaptivity in this context that I gave you is quite general. You can use it for any problems you want essentially, as long as you have a sound single level analysis, if I may say, okay. The message is that it is better in the context of multilevel to use an analysis that will go towards the L2 error, if possible, right, drive the analysis of the construction of the levels according to that. Of course, you need the bias error to control the number of levels you're going to use, but at least the distribution of the shape of the meshes should be done according to the strong error. And actually that can be proved, actually, that is optimal to do so. You just write a problem that is related to these, but instead of writing this part here, you just write the error density times the, I mean, the integers of the error density times the delta t's, you just write the whole Lagrangian relaxation, which now will involve finding functions, functions delta tL, say, from L equals 0 to L capital, and then you will see that delta tL is directly expressed in terms of the row strong, not the row weak. The only, the last level is connected to the row weak, but all the rest are driven by the row strong, and therefore that motivates directly to the use of the strong for the activity that we found later, sorry, but, I mean, the first works we did was with the weak activity. So how do we do with the time? We split in between or? As you like. Okay, I can just draw the conclusions and it would take five minutes. Let me put the conclusions and then, okay. So as conclusions, I told you again that this way of thinking for multi-hebel and adaptivity is quite general, you can use it anywhere you want. You have to be careful, though, that when you replace the statistical constraint, which means that the statistical error has to be less than tall with a prescribed probability, computationally it's much more effective to write that into a variance constraint, and this can be done as long as you have a central limit type of approximation. The central limit does not apply directly in the multi-hebel-vontical setting because you have some of independent but not identically distributed random variables, so you have to use resort to some linear or failure type of result, and you have to verify again these four each of the problems that you work with. So you have to be careful with that. Okay, so let's leave the puzzle random measures for the next one. Thank you. Five minutes? Yes, sir. Five minutes. Thank you. So we are going to discuss another idea for adaptivity in this second context. The second context involves paths that are also random and depend on time, but as opposed to the stochastic differential equations that we were given before that were etotype, these paths actually can be solved exactly. So these paths are given by piecewise constant functions, and the only random thing about them is essentially the time where they jump, and which kind of jumps we have them. These processes are very popular, they are macro processes of course, and they are used in Xenium applications, but let's go a little bit more into them. This is a logarithmic scale, this is time, and here in blue, red, and green are three components of this vector that is changing according to time. The numbers that are shown here involving are the number of species for each of these three components, so you can think of three populations randomly changing inside of a system according to maybe biological reactions or chemistry reactions for instance. Even though these green paths do not show it, they are as discontinuous as the blue ones essentially. It's just that you have logarithmic scale and you may be fooled by the logarithmic scale, that's it. So that's a motivational example, this is a gene transcription and translation you have here in this system, the component G that may react given G itself an M, M that actually may react again and give M and P, and the P that actually goes into D when it meets itself. So you see that there are sort of reaction channels, essentially when this thing happens, you know deterministically what is going to happen. The only problem is that the time where this occurs is not going to be prescribed, a theory is going to be random. And let's think of a simple example, let's think of water say in a container and water can be created from hydrogen and oxygen and it could also dissociate into hydrogen and oxygen, so this is a reversible equation. But the way the speed of these reactions actually happen depends of course on the temperature, but if you fix the temperature out of this discussion, they depend for example on how much hydrogen and oxygen is there because if there is no hydrogen, oxygen will never meet hydrogen to create water, right? So the likelihood of these events to happen depends on how dense these components are in the system. So and this can be mathematically described. So what we want to do again will be to compute some given output quantity of interest at a given final time up to a certain tolerance again with high probability. So how would describe the evolution of such a system? I told you before this was an X that lives on a lattice essentially because it's counting. So the components of X essentially describe different species in the system. So we count how much of the species one we have, how much of the species two we have and so forth, okay? So X is taking values on these lattice, right? And it was according to time and of course omega. Now we assume that we only have finitely many possible reaction channels like in the case of water for example, we may have three components, one is how much hydrogen, how much oxygen, how much water, right? And we can only have water dissociating into hydrogen and oxygen and the other two meeting to get water. So when this oxygen and hydrogen meet, right, we assume that they can only create water. So there's nothing stochastic about the fact when we know that they meet. The question is when they meet and this is the stochastic part. Okay, so when this thing happens, right, the J reaction type of in the system happens, the only thing that we have is a change from X that was the state of the system into X plus nu J. In the case of this water for example, when oxygen and hydrogen meet to create water, you have a plus one in the water and the corresponding minuses in the oxygen and the hydrogen. Now, when does this happen? Is given by the so-called intensity functions for the jumps or in this context also called propensity functions. These are non-negative functions, A, J, that are evaluated on the state X and give values on R plus. So that the probability of having a reaction of type J in an infinitesimal interval that goes from T to T plus DT, given that you have the state X is essentially A, J of X times DT. So for instance, if you have no hydrogen in the system, you know that new water is not going to be created. So this A, J corresponding to that, that reaction will be zero in that interval until that situation changes. So again, we want to compute the expected value of a given quantity for the given of several G and the applications are everywhere. So I'm not going to stop on that. So we have plenty of references and as well this is not complete, but we just wanted to show a little bit the efforts on this multilevel version. So what is interesting about these systems? First thing as opposed to the previous discussion that we have is that exact algorithms are available. What do you mean by exact algorithms in this context? You remember when we sample a la Monte Carlo to approximate quantities of expected type, we always have statistical error unless the problem is trivial and has no variance and then we are not discussing anything stochastic. So exact in this context means pathways exact, essentially that with finite work you can simulate a process that has the correct law. So there is no bias error. You can do unbiased simulation with all one cost, which means that already you have the beautiful regime that we discussed before that was total to minus 2. The cost is total to minus 2, the number of samples in Monte Carlo and then all one to solve for each path in the mean. So the whole discussion here is about constants. So you do variance reduction in order to reduce the constant in front of the total to minus 2. So that's the whole point in this context. As opposed to the other case where we actually can push the complexity from total to minus 2 log square times total to minus 2 into total to minus 2. Here it goes from total to minus 2 into total to minus 2. And still you can get implicit gain. So constants matter, that's what I'm trying to tell you. So then if you have an exact algorithm, why are you going to introduce an inexact one? Well, because the two things, because the statistical error is already there. So having a bias that is zero in the presence of a statistical error may not be optimal because the cost of simulating such exact path may be a lot. So introducing a path that is a little bit cheaper and has some bias and the bias is smaller than the statistical error will not change things dramatically from the error point of view, but it will do it from the work. So it is actually a good idea to actually introduce bias in a bias simulation. Right? Good. So what is the problem then here is that essentially the cost comes from all these events that are happening in the system. So if you actually have a lot of molecules in a big reactor, you can think that there are a lot of reactions happening. And if you have to track all of them to do this exact simulation, the cost is going to be horrible. Okay? So, well, essentially these interarrival times between jumps are exponentially distributed according to, with intensity, there is a sum of all the AJs in the system. And this can be a large number, depending on the system. Okay? So as opposed to that, then one would like to do approximate time stepping. Time stepping that is cheaper but introduces a bias. This has been done in, I mean, pretty much from the tau-lib, there is a forward-oiler method to more refined versions of it. And you have drawbacks, of course. First is you're going to introduce a time discretization error that you didn't have. Right? And two, you may have non-intuitive, at least qualitatively wrong type of behavior. Essentially, if you go back to the discussion that we had at the beginning, our process is actually counting. Counting molecules of different types in a system. These are not going to be negative, right? I mean, if you're doing the proper thing, the true trajectories will never cross zero. But if you do these sort of approximate algorithms with a non-zero probability, you may cross these things. And you may want to avoid these events because they have an issue on the accuracy. The global accuracy that you can get. So there are ways to deal with these. For example, the pre-lib is just to adjust adaptively the time step, you know, the control, the one-step exit probability. But we would like to always relate things with respect to the goal, which is to relate the old type of errors that we have in this approximation with respect to the goal, which is the expected value of g of xt capital. So what we did in this line of adaptivity was to first create a single-level hybrid algorithm that at each time step, adaptively switches between the exact path simulation and this tau-lib sort of oiler in order to do what? To minimize the computational work. So why is this somehow needed? When you're close to this zero level that you don't want to cross somehow, the tau-lib, which is the approximate forward-lowered type of discretization, does have a cost that is not being reduced somehow. So you're not gaining too much from using it. But the exact simulation essentially has a cost that I told you was essentially going up with the number of particles in the system. If you're going to get closer and closer to zero, that cost actually gets smaller. So close to boundaries in this respect, the exact simulation is advantageous because it doesn't really use bias and it's cost to do because the inter-rival times are relatively large at least. So this can be quantified and a rule that switches between these two things can be done. So again, this is adaptivity, but now not adaptivity just with respect to time, but between methods. So you're combining two methods. One that has sort of infinite accuracy or infinite order with respect to time discretization and the other that has order one, which is essentially coming from Euler. So that was one thing. The other thing is that this idea can be generalized into multilevel Monte Carlo algorithms for our control again. And not only that, it's like you can do this according to the different channels. So I told you before that when you introduce this, there are different J channels. And with respect to these Js, what matters really is the propensity. If the propensity is sort of large, you're taking a lot of events over a unit time step. If the propensity is slow, these things fire not often. So if you can actually split among the reactions, do a time splitting method and essentially treat some of the reactions with the approximate method and some of the reactions with the exact method. This can be done as well, even in the context of multilevel. And again, all this, what is going on here is that you're introducing biases and approximations that have much, much less cost when compared to the unbiased algorithm. This is why this gets useful. So how do we go about it? To motivate the different methods, the whole thing can be written through this course random time representation as accounting algorithm. Essentially the value, the number of species that you will have at time t, essentially is the number of species you have at time zero, right? Changed by what? Changed by the number of times the Js reaction occurred. And the number of times the Js reaction occurred is described here by a Poisson process whose intensity is non-trivial, it's just the integral from zero to t of A J, the propensity of the reaction J. This is again a non-negative function and this is just a random time change on a given Poisson process. So this is a natural number that goes from zero to whatever and this is just the vector that introduces the changes according to that reaction. Is that clear? So it's just bookkeeping, you see? And if this A would not depend on X, this would be trivial to simulate, right? You will just change accordingly the time in a unit of Poisson process and know exactly when the jumps are going to occur in your system and that's it. So to simplify this, when you want to do approximations, you just do a donor method, essentially you freeze the value of A at the left point and approximate this increment by a Poisson random variable with intensity A J times delta t, if you wish. And the first time I wrote it, they put tau here and this is why they call it the tau-leap method. Okay, so first problem that you see here from this approximation point of view is to find tau so that you're not going to cross these important boundaries with a given probability. Delta is given. So this tau will depend on X and delta and this will prevent you from jumping according to this inequality. What we did actually is to do some large deviation, turnoff type of bound to actually come with a computable formula that one can use to ensure this inequality. And then from this delta t, we can actually discuss whether it is worth taking the approximate step or the exact one. So and then one step farther is to go into a global probability of going out because remember this discussion that I just made was based on one step, essentially. If you're here, how you have to take delta t, so the probability of exiting in this time is below a certain value. Now what you really want is the probability of the whole path going out or not going out. So you have to link the one step probability to the whole path probabilities and then link these quantities into the effect of the global error that we really want to compute and approximate. So this is the first step that is good because I mean in principle the probability of going out for a path, if things are going to be done correctly, it's a rare event. So you want to, instead of sampling that rare event, you want to compute quantities that are sort of easy to deal with. And this approach does it actually. Okay, and like I told you before, then one can write a multilevel Monte Carlo based on these type of paths, right? Find the m's and find the l. So you control the exit type of error, essentially that is related to the indicator function of the path going out. That will produce an O1 error. So essentially that's why you have to compute to control this via the exit probability. Then the weak error that is more or less what we have discussed before with the diffusions and then the statistical error that doesn't change from the discussion that I made before. So the issue of this is here is that you have to do things correctly to simulate couple pairs that have the right law and are correlated in the right way. And of course being able to accurately estimate the different quantities in the multilevel construction. This is particularly nonstandard here from respect to the diffusions that I discussed before for the following reason. The paths, right, they're taking values on a lattice. So when things are going to converge, right, I mean as you discretize the L, the delta t or increase the L if you wish, the paths are going to be closer and closer. But you know, in diffusions closer and closer may mean you're in a continuous space, so they may mean actually getting closer and closer. But here you have a lattice. So what does it mean to get close in a lattice? You're far or you're there? So actually closer in this setting, it means that the probability of going far is going to zero, right, and the probability of being equal is getting larger and larger. What's the problem with that? That the problem is that this creates immediately high kurtosis. And high kurtosis is bad news in the multilevel Monte Carlo context because it means the estimate of variance is difficult. So one thing that appears here naturally as opposed to the diffusion case that I discussed before is you have to work a little bit farther to estimate the right variances in your system. Remember the m's are driven by the variances. So you need good estimates of these quantities. Otherwise your estimates on the m's are wrong. And if the estimates of the m's are wrong, bye-bye Monte Carlo, right? So that's important from the practical point of perspective. I think it's in delta and the index and the... Little n, yes? It's in delta that you use to balance the probability of s. Yes, everything goes there. So you know that. Everything goes into there. I mean, each level in principle here contains a delta Tl but contains also a delta L. Because remember, there was a... There was this probability of going out, right? So you may prescribe this delta L as a function of the L and also a delta L maximum per level. So this thing goes as before like delta T0 divided by 2 to the L, for example. But this doesn't go that fast. This is chosen in a different way. So these two things combine to give the method. Is that clear? So these two are varying. All right? Thank you. Good. Thank you. Okay. So this is an example, the example that I showed at the beginning. And here you see that from the exact method, essentially, which is SSA, the error bound versus the tolerance essentially goes like total to minus 2. So essentially, this is the complexity that we have. And if you use this hybrid multilevel, you can still keep the same complexity and what you gain is a constant but the constant can be substantial. All right? In this case, it's 10 times faster. And this is actually the usual plot that you would see in our papers. It's the control of the weak error. So what the user gives in this context is a tolerance but he also prescribes the confidence level. So in this case, the confidence level was prescribed to be 5%. So when you do your simulations, essentially, you expect to see around 5% of them going above the prescribed tolerance. So you see here, the tolerance is getting lower and the corresponding sample error for the different runs of the material, Monte Carlo, distributed like this and then essentially these numbers here are the percentage of those above the bound. Okay? So you change the quantity here, these exceedance values should change as well. All right? So you can actually have, depending on the problem, different type of gains. I'm not going to discuss this too much but you can go above 10, of course. You can get 10 to the 4, things like that. But the plots look like this always. You're improving the constant. You're not improving the rate. The rate was optimal to begin with. Right? So it is actually a pain to introduce a bias approximation with respect to the unbiased method. Okay? And here's an example of our reaction splitting method where you actually do this differently in the different reactions to get even farther. Okay? So this is, I think that's it. So summarizing here, this is a different type of idea for multilevel adaptive, right? You're switching between methods. Okay? You can also do the delta T adaptive as well but here you're switching between methods. And you're going from an unbiased method into a bias method to improve the constant in the complexity. All right? Any questions? Okay. So now we discuss discretization from another point of view. We're going back to uniform discretization but now we ask a different question. We say we have been driving our discretizations with this type of geometric type of refinements, right? I mean you get an initial discretization parameter, then you divide it by a given number, then you divide it again by the same number, and so forth. Okay? This is natural within, again, numerical discretization, right? You try with a mesh, then you start to divide it. Try with another mesh. It's not optimal. It's not optimal. That's a natural question to ask, right? I mean, should we choose the HLs or the delta TLs like this or should we not? Answer? No. It's not optimal. Okay? So with some assumptions you can show that this is not optimal. This is the work we did in this optimization of hierarchies for Matias Monte Carlo. And this other work that is here, this continuation of Matias Monte Carlo addresses this problem, right? That these variances that you need actually to drive and get the optimal ms are not known a priori and somehow have to be found on the fly, right? So you need to create an algorithm that not only is adapting the ms but is trying to find these parameters that are needed to find ms, okay? So that's what the continuation is doing. And this is implemented in the lecture in this code, MIMC-LIVE, this library, MIMC-LIVE. So let's go to the assumptions then to focus. So what do we have here? We have a Matias Monte Carlo method as before, but now the levels, GLs, are going to be driven by these HL parameters of discretization and we do not say that they are geometrically distributed a priori, okay? So we're going to let them vary. So we're going to take L capital plus one of those and we are going to assume that essentially the weak rate is algebraically convergent. The strong rate again is algebraically convergent. And the work, for example, is algebraically divergent. That's in the linear PDE that we have before. And then what we have to find are not just the ms but also the Hs, okay? Is that clear? So as you know, the constants and the Qs here will depend on the problem and the solver, as we discussed before. And these are just a couple of examples, right? In Euler-Maruyama you have Q1, Q2 equals one. And if you have the other PDE example, we already discussed that the Q1 was two and the Q2 was four. Gamma again is defined by the solver and it's obvious in Euler-Maruyama, it's just one. So I don't stop into that. I just try to show you that these things can be found. And if you have no theory at all, right, if you're desperate, you can just try to fit something to the convergence as you see, right? And or the explosion in the solver. So you can actually measure things as well. So two useful quantities here is the quotient between the strong convergence and the explosion of the work that is chi and the eta is the quotient between the weak convergence rate and the explosion of the work. Of course, the best case for multivariate Monte Carlo is when chi is larger than one, the limit case is chi is equal to one. And the sort of algebraic improvement with respect to Monte Carlo is chi less than one. All right? And again, you see that if you change the solver, you may change from chi greater than one into chi less than one as I described here. In the PDE example, if you use an iterative solver, you get 1.34 and if you go into the direct solver, you get 0.89 already. So you can influence the convergence rate of multivariate just by changing the solver like that. Okay. So the optimization, how does it look? It looks like this work that was the work that you discussed before, subject to this constraint. This is a bias constraint less than this fraction of the tolerance. And this is the statistical constraint, which means that the total variance, right, is less than theta times tau divided by c alpha. And c alpha, again, is the confidence constant that is coming from the user, essentially, when he imposes the confidence level alpha. Okay. So the good news is that you can solve this problem. And of course, you can even verify the normality, how far you are or how close you are when you do actual computations. And here's a result somehow for chi equals one. You actually have, is a case where you actually have geometric discretizations. But I have to warn you, these geometric discretizations, you see, h a is somehow beta to the capital minus l. So it's a subdivision that we discussed before. This doesn't depend on the levels, but it depends on the total. So it's not your usual geometric discretization that you would think of. Okay. And the beta, actually, is not two or three. I mean, it's something that may even have an influence from total. Okay. So the spacing between the levels is not trivial either. All right. This is one more thing that is interesting here. The optimal choice of the splitting parameter, which is the fraction that you devote to the bias versus the fraction that you devote to the statistical error, actually goes to zero, sorry, goes to, oh, come on, goes to, I mean, the fraction of the bias goes to zero as a total goes to zero. Yes, that's right. So essentially, if you see here, theta is the fraction that you give to the statistical error and one minus theta to the bias. So theta essentially goes to one, why, because theta is fixed and L goes to infinity. So this theta is just going to one. All right. So asymptotically, you're putting most of your error, you know, your error allowance to the statistical error, not to the bias. And of course, the level behaves like log total and with this constant. But for non-Ka equals one, things become less simple. You can still find the relations on the theta and you can find bounds on the number of levels with respect to the different parameters as this shows. Okay. Nothing really changes from the classical results in terms of exponents, but you can identify how the input constants affect the actual complexity in the actual constants. So you can play tricks to reduce these constants knowing that they exist and the shape they have. Okay. So the shape they have in these constants in the different cases are like this. So you can actually see that in some of them, for example, the variance of the level zero plays a role, but in some others they don't. So here in this setting, if you're here, you may get an advantage by doing a special trick for level zero. So whereas in the other cases, you may not. Okay. So like I said, the exponents are the same as you have in the classical theory, but now you identify the constants as well. And you identify how you have to space the levels too. And that is more or less the discussion that I have here. These are not going to be necessarily nested because you have this beta that behaves in a strange way. And the H0, also the first mesh, also behaves in an interesting fashion too. Here's a discussion on the continuation of multi-level Monte Carlo. And again, the plots always look like I discussed before. You have a direct solver here and an interactive solver. You have the prescribed error. And here's it versus the tolerance. And you're checking essentially that when you run over and over and over the multi-level approximation, you get a fraction of those that go above in accordance to the more or less the confidence level that you prescribed. Now coming back to this question of finding the variances and other relevant parameters to do optimization, in the continuation of multi-level Monte Carlo, what we do is to solve a sequence of problems until you hit the problem that you want. So essentially, what happens? We know that in the best possible case, the explosion of the work in multi-level Monte Carlo is told to be minus 2, right? Best possible case. It could be even worse. And you can have a long talk here. What that mean? The worst is like this. So if you actually solve now a sequence of problems where you first begin with told zero you'll say, and then you go times told zero divided by a certain fraction, right? Call it F, and then keep going, right? Until you hit the actual tolerance you want, you can sum these works here, and you can get a choice on this fraction so that the actual sum of these sort of preliminary works does not create a huge overhead on the last one that you want to solve. That's the first idea. And then what happens here? That when you start solving these sort of axiomary problems, you start learning the constants that you need, okay? So somehow by creating these fractions, you do not overshadow the actual work that you want to do optimally, but you start sort of uncovering a veil over the variances that you need and the other things that you may need to find the optimal elements. That's essentially what is behind the continuation of Monte Carlo, sort of a motto with respect to a tall type of idea. Not really complicated, but it has to work quite nicely. Of course, on top of that, there are tricks to actually improve a little bit the variance estimates on the deeper levels, right? What happens statistically is that in the course levels you have a lot of samples, just by construction, right? And Monte Carlo produces a lot of samples in course levels and very few in the deep ones. But the deep ones are the really costly ones. So if you make a mistake in the variance there, you make a mistake on the M and things go really rotten. So to get a better hold on the deep levels, you actually do some kind of Bayesian type of estimate where you use as a prior model an extrapolation that is coming out of the course levels and you enrich it with a few samples that you have from the deep levels. That's what is behind too. So it's not just one idea, but I mean, I'm just telling you the couple that I consider most relevant. So the good news again is that if you do all this sort of little learning process, the complexity that you get is not much of an overhead with respect to the optimal line here. That will mean that for this sort of pink magenta type of a curve, you are just writing the cost that you would incur if you would know all these nice constants beforehand with respect to finding them on the fly. So it's not too much that you have to pay, really. So that's with respect to the optimal splitting theta and it's not really much relevant. So what is the second part of conclusions? We didn't discuss much of normality, but we have a result using the linear failure in this context. Better splitting between Bayesian statistical error can improve the work a little bit and this continuation of Monte Carlo tries to provide a more stable algorithm that allows you to find the constant in between. Geometric hierarchies are nearly optimal, but they are not the geometric hierarchies that you are used to work with. They come from this optimization. They have specific values that may depend on toll as well. So for ETOs, I included in the slide some assumptions and some results with respect to the densities that I don't have time to comment on, but I'm going to share the slides so if you have questions you can come to me. Again I emphasize here that these files actually depend on derivatives of the final solution with respect to initial conditions. So this is what really makes the method non-adapted to the natural filtration. So let's say this is multilevel in that context. This is what I meant by strong convergence and the accuracy result. This is what I meant by the work optimality result and I included a little bit more in terms of complexity. But OK, CLT2. So last thing before I run out of time, a little bit of an activity with respect to M and I just touch upon the classical expected value computation using just Monte Carlo just to close it up. So what's the point here? You want to compute the expected value of X using IID sample averages and you want to prescribe again this delta which is the confidence level and you want to prescribe also the toll. So what is the point here? The point is that usually you do this approximation based on the central limit theorem. That rule, that will be a rule of the following type. Essentially you are going to prescribe toll and you are going to prescribe delta and then a number of samples M will come out. It may not be so reliable when you have X that has low decay and density, it's essentially heavy tails. So for that reason we provided a rule that is a little bit more robust and essentially combines the central limit approximation with other type of bounds. Let me walk you through a little bit before I finish. I have 10 minutes, right? Yes, okay. So this is going to be a little bit intense then. So this is a general stopping rule, right? Based on the samples you are going to estimate the probability of failure and based on that we are going to accept or not. So the probability of failure actually is the probability of your approximation not being close enough to toll with the prescribed probability. So this event of failure having a value larger than delta essentially. So again as we described with the central limit theorem before the use of variance constraint you can actually, if you believe that your sample average, your sample variance is actually accurate enough you put a rule of the central limit type substituting sigma here which is a true standard deviation by its approximation and if this actually is accepted then you stop and you give your approximation. These may have issues, okay? So consider the following Pareto distribution. The density is decaying like X to the minus 4.1, right? For X larger than 1. And you see here the behavior of this quotient, this quotient is the probability of the actual failure versus delta the one that you wanted to impose. So this quotient is larger than 1, this is bad news, this quotient is less than 1 or equal than 1, it's okay. So you see here tolerance versus delta, right? So here it's like having more and more confidence in your computation and here is more and more accuracy as you go to the left, right? So as you let tolerance go to zero essentially for delta fixed things are going to be fine but in the regime where delta is say a few percent and sorry, when total is few percent and delta is sort of stringent and that's natural for engineering computations this is not going to be nice, it's going to be actually in the red and it's not going to give you the prescribed tolerance. So there are many results in the literature where you fix the delta and you let total go to zero. What I claim is that essentially those results overlook this regime which is important in the applications. Okay, so the whole point here is to penalize this fact because we really don't believe sometimes that we have such an accurate estimation of the variance to uncover this what we do essentially is to use very same type of bounds that are essentially sort of better control rates on the approximation in the central limit theorem and combining those with the sort of H-worth type of expansion which is something that you would use when you have full confidence that the tails are going fast essentially to zero then you get sort of an improved rule and when you go above this rule you actually see in this example that you can achieve the right prescribed failure rate with almost no change, almost no change you see in the computational cost and the cost here is just determinants of samples because we are working with just Monte Carlo in this context. And this can be done also for another kind of not so nicely behaving random variable and you see again that for small delta and relatively few percent total this is a bad case but you can fix it using this algorithm and without incurring too much in an increase of the work. And finally this is an example where we actually expect the usual approach of central limit theorem to work nicely and it does actually the work here is such that the failure probability is below the delta but this algorithm also achieved that with that an increase, a sensible increase in the work. So in the nice cases this behaves well, in the worst cases or the bad sort of cases this will have better control on what is going on. Okay so that's a summary and I think I get to my end that's remarkable actually. So thank you. There any questions or comments? Yes, Amar. So I would be interested to know if we could apply this criterion for example in the case where X is related to Rar event. So that is a Bernoulli random variable typically which is 1 and 0 but 1 with very small probability. What? That's not the case we looked into. I think that needs a measure change. I don't think one can apply it directly here. Just using some simple averaging of your output. So our robust is the criterion to... I mean what I'm saying is that it can be used but provided the first you do a measure change. If you don't do a measure change all the estimates are going to be wrong, right? You mean wrong inaccurate in fact? Absolutely. The whole thing is completely reliable, right? So Rar event asks for a measure change to start with and then if the measure change is sound, the important sampling is sound then you can use things like this. But I wouldn't recommend using this in a Rar event without looking into an important sampling algorithm first. It can be combined with the same... That's what I'm saying. But I wouldn't recommend it to use it without the important sampling step first. I have a second question perhaps. Because again that regime is different, right? I mean that regime is like for example 0, 1 with a very tiny probability of achieving 1. I mean here it's like unbounded, right? You can get unbounded variables with a certain probability and this is the case, not that fast. So you can achieve them with the probability that it's not so negligible and this is what gives you trouble in the estimation of the basis. It's not the same regime as the... Have you looked in the direction of concentration of measure inequalities like we know... Well, that's also not applicable here in the sense that for concentration of measures at least if you use classical large deviations you need exponential moments. These guys do not have any, right? I mean we use them in the Cherno for the Taulib but here you don't have an exponential moment. So if you have a concentration it's too little. A short comment. Maybe the applicability of your method is connected with a kind of ratio between the expectation and the variance. Because the case of Fourier events is that 1 over square root of p goes to infinity and p goes to zero, something like that. The relative value of the expectation... Yes, yes, but I mean that's also what defines the regime of red events, right? Yeah. But in this case again I emphasize that is the tail that is creating the issue and we create sort of something that works more or less uniformly with respect to that and does not address at all the question of the red event. The red event should be addressed in a complementary way. At least it's my humble opinion. Other questions, yes? Have you looked or is there any chance to extend this framework to the case of PySYS adaptivity? I'm thinking of SPDE where you would like to have a discretization that has some local adaptation in space. But this is not about this last part, it's the previous part. Yes. Well, I'm kidding. Okay, okay, okay. This is about the... Again, the adaptivity would respect to... Pathwise. Pathwise, okay. So what we discussed was pathwise but in time. Yeah, but the time step is fixed. No. Independent of the... No, no, no, no, no. The time step is delta T of T and omega. Delta T is fully stochastic. So it cannot have two features that come with omega, pathways. So yeah, that's actually the most general but you only use it if you're obliged to. That's my advice. That if you can use, for example, the terministic type of activity, that is cheaper. But you need some prior information about the problem that you're solving. Right? Okay, thank you. Did you have any photo questions? If not, so let's thank again the speaker. Thank you. Thank you. Thank you. Thank you. Thank you. Thank you.
the time allows, I will also discuss optimal use of your favorite method. In this context, you corresponding complexity, which is the computational work that is needed to achieve a for a prescribed confidence level, okay? And to by this, the statistical error governed by that. Giles. Let's take a constant beta for differences between consecutive approximations. is the difference between two random variables. here, this is saying at the beginning I only had essentially, thanks to the construction, based on is essentially saying that in the deepest level of the subdivision multiplicative constant beta. right? So essentially both the work and the constant here. So the total result here is L on introduction with the uniform methods. Now the ideas to introduce non-uniform discretizations or time steps such as t0 is 0 and the final time least to show results on the algorithm in this known a priori, right? Because this is the number refinement, you need to produce essentially the Bridges, right? That is the correct way to do asymptotically you have normality, asymptotically discretization here is the accuracy constraint, asymptotic normality, you also have asymptotic between uniform Monte Carlo and adaptive Monte message is that it is better in the context of the analysis of the construction of the levels error to control the number of levels you're algorithms are available. What do you mean by approximate time stepping. Time stepping that is again, this is adaptivity, but now not adaptivity Monte Carlo algorithms for our control again. And uniform discretization but now we ask a different optimal choice of the splitting parameter, which the first mesh, also behaves in an interesting Bayesian type of estimate where you use as a of the course levels and you enrich it with a few number of samples M will come out. It may not be have such an accurate estimation of the variance