WEBVTT 00:00.000 --> 00:11.920 Starting the next session five minutes early. So to give our next speaker, Lisa, this 00:11.920 --> 00:19.440 is our own Piscuit policy, no, this Piscuit, hold on, I'm sorry, Piscuit cobble, it's 00:19.440 --> 00:24.440 a tool, okay, apologies. Just to give him five more minutes, because it's a really 00:24.440 --> 00:30.920 depth talk, which is basically velocity sampling efficiently from high dimensional distributions, 00:30.920 --> 00:41.720 I would say you. Okay, thanks a lot. Okay, so today I will try to present this topic, 00:41.720 --> 00:50.280 I will give some brief background and then I will present the software that we have. So let's 00:50.280 --> 00:57.640 start of understanding the title. So what is a distribution, and distribution in let's say 00:57.640 --> 01:03.320 two dimensional spaces, we can have a Gaussian distribution. So for example here the points are 01:04.600 --> 01:12.280 a Gaussian distribution, and we would like to study truncated distribution. So you have points on 01:12.280 --> 01:19.640 the plane that are in Gaussian, and then you have to truncate with some shape. This shape can be 01:19.720 --> 01:26.280 convex or non-convex or linear, non-linear. So here we have just a polygon. So 01:28.120 --> 01:35.080 another distribution, so this is in one dimension, we can have a function, this is a function 01:35.080 --> 01:44.360 that describes the distribution, which is we call it logoncave, and then we want to sample points, 01:44.920 --> 01:51.000 we truncate this distribution, and we want to sample points from minus one to one. And this is 01:51.000 --> 02:00.600 the histograms that we compute. So imagine this objects in with more variables, then we have a 02:00.600 --> 02:06.840 high dimensional distribution, and we want to sample from this. So what is the main, so this is 02:06.840 --> 02:11.640 this are the objects that we want to study, and the problem is that we want to sample from this 02:11.720 --> 02:17.800 distribution. So let's say that you want to compute points, if you have a uniform distribution, 02:17.800 --> 02:23.480 you have a polygon, and you want to compute points. This is a simple example, you want to compute 02:23.480 --> 02:31.560 points that are filling your polygon. So apart from being the fundamental problem in mathematics 02:31.560 --> 02:37.960 and computer science, this sampling is a building block for if you want to do Monte Carlo with 02:38.040 --> 02:46.600 a grayion and volume computation. So why is this interesting, because there are many, many applications, 02:46.600 --> 02:55.080 so for example, in bias and inference, when you want to estimate parameters that have constraints, 02:55.080 --> 03:02.200 so the constraints now is this truncation that we show, your polygon, gives your constraints, 03:02.840 --> 03:09.400 or if you want to do constraint optimization, this is a well known example is a linear programming. 03:09.400 --> 03:15.000 So you have a poly, a poly top, and you want to optimize a function over there. 03:17.640 --> 03:26.120 Another area of application is in finance. In finance, we have, when you have portfolios, 03:26.200 --> 03:36.360 or you play stocks, then you have constraints. So constraints come from some market, or constraints 03:36.360 --> 03:46.920 of constraints over the portfolios, they can be probabilities, and an arbitration is 03:46.920 --> 03:52.280 incapable of biology. So in metabolic networks, when you want to study the metabolism of 03:52.360 --> 04:02.200 organisms, then interestingly, you can model it with the high dimensional poly top. This is 04:03.720 --> 04:10.440 high dimensional example of a polygon, and then the biology is going to sample from there, 04:10.440 --> 04:17.800 from some distribution, and this gives some information about the organism that you want to 04:17.880 --> 04:23.720 study the metabolic. So for finance and biology, I will give certain examples. 04:25.000 --> 04:31.240 I hope that this will be more clear than it's next slide. So if we have a simple, 04:32.840 --> 04:36.840 object, a simple geometric object like a triangle, or, 04:38.520 --> 04:44.440 in fact, the measure we call it a simplex, or we have a hyper cube, then it's easy to sample from 04:45.000 --> 04:49.560 there. We know how to sample. And you can do it also in high dimensional. So they're 04:49.560 --> 04:56.040 formulas for this. You just evaluate the formula and you have your points. But I assume that you have 04:56.040 --> 05:06.680 something, a more strange shape. One way to sample from it is to compute the bonding box, 05:07.240 --> 05:13.560 and then sample from the bonding box, and say, okay, I will sample from the bottom box, 05:13.640 --> 05:18.280 and then see how many points fall inside my object that I want to sample. 05:18.280 --> 05:24.440 Okay, into this is fine, right? So you sample points, all this red and green points, 05:24.840 --> 05:30.040 and say, okay, how many points are the green ones? How many points are all of them? And then 05:30.040 --> 05:38.840 you can estimate the volume of your unknown object, or you can sample from it. It's the same 05:38.920 --> 05:48.760 problem. Okay, this shape, after a few dimensions, so if you draw this in 15 dimensions, 05:48.760 --> 05:54.040 you will not be able to compute the green point. You will sample from the box, 05:55.080 --> 06:01.000 and all points are red. Okay, so in high dimensions we have to do something 06:01.640 --> 06:12.280 and more clever if you want to sample from this unknown shape. Okay, and the solution to this 06:12.280 --> 06:19.480 is random groups. What's a random group? You have a point inside your convex set, 06:21.000 --> 06:26.120 the set that you want to sample from, and then you can move to a neighboring point, 06:26.120 --> 06:36.200 and you continue this procedure for many steps. Okay, so this can be analyzed as a 06:36.200 --> 06:42.040 Markov chain. I will skip all the details about analysis, and if this chain 06:42.040 --> 06:51.080 converted this or not, but just as a general idea, a random group is computing points, 06:52.040 --> 07:02.040 like walks inside your polygon. Okay, and the two examples here is, this is a two 07:02.040 --> 07:10.520 dimensional polygon, but you can do this also in non-convex shapes. Non-linear shapes, 07:10.520 --> 07:18.440 that are convex, you can also do it in non-convex shapes. Okay, so let me now give some examples 07:18.440 --> 07:28.760 of some walks that we use also in practice and these are implemented in our package. So 07:30.040 --> 07:39.080 I will not discuss about the algorithm, I will try to explain it in simple words. So assume that 07:39.160 --> 07:49.720 you learn how to play a billion. When you start playing a billion, you take, you can hit the ball 07:49.720 --> 07:59.320 with not a lot of power, so because you are learning billion. So your ball will be moving inside 07:59.640 --> 08:07.560 the big table, little by little, so you are learning how to play billion. So this is 08:09.400 --> 08:15.000 ten meters of this work, this is the first work, the most simple work, which is called ball work. 08:15.000 --> 08:23.560 So you have a point, you compute a ball with some certain radius, and then you move randomly 08:23.640 --> 08:30.280 in this ball with some probability. If the, if the new point is inside our convex body, we keep 08:30.280 --> 08:36.600 it, otherwise we throw the way. And then you continue likewise in the next step. So you are here, 08:36.600 --> 08:43.080 you compute again a new ball, you sum randomly, and you move. So this is a ball work, this one basic work. 08:44.120 --> 08:51.080 You probably know maybe you know the Metropolis casting algorithm, which is very 08:51.800 --> 08:59.160 has a long history. This is the variant of ball work with our constraints. So if you do this 08:59.160 --> 09:03.000 without the constraints, if you have a point, you compute the ball and you move there, some probability, 09:03.000 --> 09:10.680 this is a Metropolis casting algorithm. Okay, now assume that you are becoming a better player in 09:10.680 --> 09:22.440 billion, and you want to, to suit the balls with more power. Okay, now you are shooting in lines, 09:24.200 --> 09:32.120 so you can suit the ball further away in your table, you are more confident. So this reminds us of 09:32.120 --> 09:40.280 the next work, which is a hidden run. So you are in a point, inside your convex body, 09:41.160 --> 09:49.080 you compute a random line, this randomly, throw randomly from a ball here, a random direction, 09:49.640 --> 09:56.280 this defines a line, you take the segment, which is inside your bolt open, then you move on this segment. 09:57.000 --> 10:05.160 Then you, you draw another line randomly, and you move there on this segment. 10:06.120 --> 10:13.240 So this is a more advanced way to move in your bolt up, because you can go intuitively, 10:13.240 --> 10:20.600 you can imagine that you can go further away. So the work can converge faster, so in practice 10:20.600 --> 10:29.720 converges faster. Okay, and now that we can assume that you are even more advanced, you learn 10:29.800 --> 10:40.200 billion even better, so you want to do something like this. You want to throw your ball, 10:40.200 --> 10:53.000 and then also hit in the borders and do some trajectory reflections. Okay, this is called the 10:53.000 --> 11:01.320 billion work. So you have a point, you select the direction, randomly again, and you suit it there 11:01.320 --> 11:08.600 with some length, so the length is a parameter of this work. So your next point will follow the 11:08.600 --> 11:19.480 trajectories and end up in some other direction. In some other position, inside your 11:19.480 --> 11:27.320 complex body. Okay, this is called the billion work. Okay, just a small comment here that 11:27.320 --> 11:36.920 there are even more advanced works. When where you, you can hit your, you can move inside the 11:36.920 --> 11:44.520 body with some curvature. So here all the, all the trajectories are linear, but you can also 11:44.600 --> 11:57.000 follow some care of the trajectory. And the name of the work is the Hamiltonian model, 11:57.000 --> 12:02.040 the Carlo work, which is based on Hamiltonian dynamics. I will not present it here. 12:02.600 --> 12:16.280 Okay, now one question that comes in mind is how many steps should I do in order to have a point 12:16.280 --> 12:24.840 which is random. So this is an experiment with the works that I present to you. The first row 12:24.920 --> 12:30.760 is the ballwork, the second and the third is the hit and run, and the third is the billion 12:30.760 --> 12:41.480 work, which where you have also reflections. In the beginning we have our work with a very small 12:41.480 --> 12:50.040 step, you can see that the ballwork cannot conferred, but if you give them more steps for a point 12:51.000 --> 12:57.800 to convert, then your sampling is better. So here we try to sum it uniformly from a high 12:57.800 --> 13:05.320 dimensional cube, the cube is 200 dimensions. And this is just a projection in 3B. So hit and run 13:05.320 --> 13:13.560 is a bit better, and interestingly the billion work, even with one step, can convert. So even 13:13.640 --> 13:21.560 one step of a billion work and have a point that is well, you can see it as a random point. 13:24.680 --> 13:31.400 Okay, but in general, let's say that you implement this. It's not possible to draw pictures and see 13:32.680 --> 13:39.480 if our points are nicely distributed. Of course there are theoretical bounds for all these works 13:39.560 --> 13:46.920 or for most of them. That tells you, okay, at the point, if you run this random work for 13:48.200 --> 13:55.800 some bound on the number of the dimensions, the point will conferred, but as you many of you 13:55.800 --> 14:01.800 probably know, in these theoretical bounds, there are huge constants, so they are not practical. 14:02.920 --> 14:08.840 What we do in practice, we can apply some statistical tests, like the effective sum of the size, 14:08.920 --> 14:15.960 or the potential scale reduction factor. So these are test from statistics that you can use to see 14:15.960 --> 14:23.400 if your points are randomly distributed according to the distribution you want to sum. 14:24.600 --> 14:30.920 And a challenge here is that if you use the sum link as a subroutine in some other algorithm that 14:31.000 --> 14:38.440 builds on sampling, then it's not clear how you provide this statistical guarantees. 14:39.640 --> 14:46.120 So these are tools that have theoretical guarantees, those guarantees cannot apply it in practice, 14:46.120 --> 14:59.160 so in practice you do tricks. Okay, so these are these arguments are implemented in a 15:00.120 --> 15:10.040 package called the velocity. It is a C++, and there are two interfaces, one in R and one in Python. 15:11.400 --> 15:17.000 There are different vengeance, VR packets, maybe it's the most mature that it's, it's also, 15:19.400 --> 15:26.600 you can find it in cran, which is the repository for the packages, VR packages, so you can download 15:26.600 --> 15:34.360 it from there and use it. It contains algorithms for sampling, but also, I'll use for volume 15:34.360 --> 15:44.840 integration that builds on sampling, and copul estimation, so copulas are, are these 15:45.640 --> 15:53.800 joint distributions between two or more random variables. So this is how copulas two random 15:54.120 --> 16:06.840 variables look like. This is taking from, from the biologists, this metable with the 16:06.840 --> 16:15.880 biomass of the, of the body. So I'm not explaining it, just to show how copulas look like. 16:16.200 --> 16:23.720 And then it supports optimization for different body, so apart from polytops, which are linear, 16:23.720 --> 16:29.960 you can have a spectahitra, which are, both like this, that appear in semi-defin programming, 16:31.160 --> 16:39.320 and you can have also zonotops, that appear in some applications in robot motion plan. 16:39.320 --> 16:48.120 OK, and the last thing that is contained in this package is utilities for financial and biological 16:48.120 --> 17:05.640 applications, that I will briefly explain later. OK, so now I can briefly show the art in the face 17:05.720 --> 17:18.680 of, of velocity to have an idea how you can use these functions more or less. So here is the, 17:18.680 --> 17:26.520 the example from, something from the square. OK, this is the example of I, I created in this light. 17:27.800 --> 17:33.960 And this is the experiment that I told you that if you go to high dimensions, then you cannot, 17:36.600 --> 17:44.280 you, you cannot compute the semi-correctly. So this is just a, a simple script that computes, 17:47.000 --> 17:53.880 computes points on the square, and then asks if this point falls inside the, the ball, 17:53.880 --> 18:00.680 but goes to, to hire dimensions, as you can see, 18:01.320 --> 18:11.960 into dimensions. So here is the, the dimension, the volume of, the, the, the circle, 18:14.280 --> 18:20.840 the exact volume of the circle and the, the ratio between those two. So you can see that if you 18:20.840 --> 18:28.520 will increase the dimension, up to dimension, this is the, the error starts increasing, and then at some point, 18:28.600 --> 18:34.680 you will not find, in dimension 15, you cannot find the point inside your, uh, certainty. So, 18:34.680 --> 18:42.200 this, rejection sampling, rejection sampling cannot work, if you go to, let's say, more than 10 variables. 18:44.200 --> 18:53.560 OK, so let's see how the, implement the drawdown works are, uh, solving this problem. So, 18:54.520 --> 19:03.240 this is the experiment, uh, with, um, uh, sampling from this high dimensional cube that I, 19:03.240 --> 19:09.960 showing in the slides, uh, here we select, we can select, uh, the work and the work length, 19:10.680 --> 19:17.000 how many steps the algorithm is doing, and, uh, we can have, for example, 19:17.960 --> 19:27.560 uh, if the work step is small, you see that the work is not converges, uh, otherwise for a longer, 19:27.560 --> 19:35.640 uh, step, you have some convergence, and with a long step, uh, you have converges in there for the 19:35.640 --> 19:47.080 bold work, um, this and experiment, we, OK, this is running a sample for, uh, 100 dimensions, 19:47.080 --> 19:54.200 you see that it's quite, uh, fast, uh, so I sample 1000 points in 100 dimensions right now, 19:54.200 --> 19:59.720 uh, with billion work, and this is the result, uh, with only one work step. So, 19:59.880 --> 20:08.360 this converges, uh, really fast. Uh, of course, the cube is not, it's, it's a simple, uh, 20:09.640 --> 20:13.400 convex region, but, uh, it's good for, uh, the demonstration. 20:15.080 --> 20:23.400 OK, and then this is the, uh, I can run all of them, and, uh, uh, compute some statistical tests, 20:23.480 --> 20:30.360 so this is the ESS test. It tells you more or less how many points are, um, 20:31.480 --> 20:38.760 how many points are statistical independent. So, let's say here, we compute 1000 points in 100 20:38.760 --> 20:46.200 dimensions, uh, with hidden run, uh, only three or 19 are independent. Let's say that random 20:46.280 --> 20:52.280 is distributed nicely. Uh, in bulk work is even worse, and in billiard work, as we also saw, 20:52.840 --> 20:57.000 uh, before you have a very big number. So more than half of the points are, 20:57.000 --> 21:05.640 nicely distributed. So this gives you a quantitative, uh, test of the pictures that we show before. 21:06.280 --> 21:12.920 So we can use it in practice. So we can implement, uh, an algorithm and, uh, computer 21:13.080 --> 21:18.040 that's a, this is called test, and if this statistical test is, uh, good, we, we continue. So we can 21:18.040 --> 21:28.600 actually use it. It's not just like a plot. Uh, OK, and then, uh, we can, I, I'm not describing the 21:28.600 --> 21:37.800 algorithm, but using sample, you can, uh, you can implement, uh, algorithms for, uh, volume computation. 21:38.280 --> 21:45.480 So let's say that we want to compute the volume of, um, of a cube in four dimensions. 21:46.440 --> 21:53.640 We know that this will be 16. Uh, so this randomize algorithms. If I land many times, I get different 21:53.640 --> 21:59.960 values around around 16. So I can do, uh, statistics there that can increase the, um, 22:00.040 --> 22:10.440 varied bound. Um, but I only want to show this. Uh, so here I can, I can, I can, I can compute statistics 22:10.440 --> 22:19.240 on the, I can compute many volumes and then do some, uh, mean and median, uh, uh, 22:19.240 --> 22:27.160 computations to have some statistics. Uh, I only want to show one example here that shows the, 22:27.240 --> 22:33.240 the power of these methods. Um, here is a well known, uh, a portal from that, 22:33.240 --> 22:39.880 combinatorics. This is very difficult to compute the volume. Um, the 10th version of this 22:39.880 --> 22:49.320 polytope is, uh, in 81 dimensions and the volume is this number, this fraction. This is computed 22:49.400 --> 22:59.960 by a massive parallel computation after, um, I think weeks. Um, and, uh, this is the exact volume. 22:59.960 --> 23:07.320 So we can run our software on this polytope to have an estimation. And this, uh, only takes, 23:09.640 --> 23:16.280 some, uh, seconds. And we have something that hits the order of magnitude. Of course, we 23:16.360 --> 23:26.840 kind of expect to have, uh, the same number. Uh, but, uh, we can, uh, we can use it for, uh, up to 1000 23:26.840 --> 23:34.040 dimensions and compute the volume of these polytopes in 1000 dimensions. So, this, so, uh, 23:34.040 --> 23:39.160 showcase the, the scaling of, uh, these methods. And, uh, you can also use it in the 23:39.160 --> 23:45.320 regression. So you have a polynomial and you want to integrate it over, uh, a convict set. 23:46.200 --> 23:52.920 Uh, this can be either, this can be done by this packet screw back to, but since this is deterministic, 23:52.920 --> 23:59.720 it cannot go to high dimensions. Uh, I have an example here, uh, with the list that can compute, 24:00.040 --> 24:09.480 uh, a good estimation, uh, a good estimation of, uh, this is Degral. And there is an interesting 24:09.480 --> 24:18.280 Python, uh, packets by, uh, European space agency that computes, uh, this kind of integrals, 24:18.840 --> 24:31.720 uh, with the GPUs. So I didn't test, uh, the comparison, but, uh, they, they try to, to, to overcome 24:31.720 --> 24:39.480 the, the problem of, uh, uh, uh, uh, high dimensions with the GPUs. So this also, 24:39.800 --> 24:46.600 any interesting, uh, direction. Okay, let me go back to slides to slides. 24:47.560 --> 24:56.520 And just so the two applications that, uh, uh, we can have, uh, with sampling. So first, uh, 24:56.520 --> 25:03.080 finance allocation. So the set of portfolios, uh, with this investment on a set of stocks, 25:04.040 --> 25:10.760 uh, can be represented as a simplex. So, uh, if you have two stocks, then, uh, you have a triangle. 25:11.400 --> 25:18.280 And all the possible portfolios are points in this triangle. Okay. So you can also have extra 25:18.920 --> 25:25.880 constraints, uh, so you take this triangle and you cut it by new constraints like, okay, I want to play 25:25.960 --> 25:32.600 stocks, only, uh, U.S., or whatever. And you create some polygon. So imagine this, if you 25:32.600 --> 25:39.240 don't have only two stocks, if you have thousands of stocks, then this is a huge, uh, uh, political. 25:39.240 --> 25:46.440 And you want to study this, uh, we started this particular scenario where portfolios have the 25:46.440 --> 25:52.440 same volatility. So volatility is the degrees of a variation of, uh, of the prize. So 25:53.400 --> 26:01.000 this, this means that all the stocks with the same volatility, they are, are, they lay on the surface of 26:01.000 --> 26:09.000 a ellipseoid. So into these, uh, is an ellipse. Uh, so you want an intersection of these ellipse 26:09.000 --> 26:17.240 with a triangle. Uh, as you can imagine in, even in three dimensions, this shape is crazy. It's non-convex. 26:18.120 --> 26:25.240 Uh, and, uh, we, we implemented some, uh, algorithms based of these sampling 26:25.240 --> 26:32.920 angles that I explained you before with random walks, uh, that started, uh, uh, the anomaly detection, 26:33.480 --> 26:40.520 in those stock markets. So the talent is to, to sum up from this, uh, objects that are not 26:40.600 --> 26:48.200 complex and in high dimensions. Okay. And, uh, an application in structure biology, 26:48.200 --> 26:56.520 interestingly, again, uh, if you go to, to study with a volcanic networks, uh, again, you have, uh, 26:57.720 --> 27:06.200 a polytop or into the, you have a polygon. And the points, uh, there are, uh, um, 27:06.280 --> 27:12.440 specific cases of, uh, uh, of the flow in your network. Let's say, 27:12.440 --> 27:19.080 kind of adding that there is a flow in the metabolic network, because, uh, they, they are, 27:19.080 --> 27:29.400 connections between, uh, metabolites. Uh, so, uh, you, you want, you, if, if the biologist 27:29.480 --> 27:37.640 is going to study the balance state, uh, the state state of, uh, this, um, uh, this system. 27:37.640 --> 27:45.080 And this is a commercial job. So you can use exactly the same, uh, algorithms there, uh, 27:45.080 --> 27:53.960 and, uh, you, you can study, uh, for example, uh, uh, uh, uh, uh, you, you can start the 27:54.040 --> 28:01.480 some reactions and, uh, compute histograms like this, uh, by sampling. So, this is an interesting 28:01.480 --> 28:07.480 challenge because those polytop are, uh, since the networks, those metabolic networks can be 28:07.480 --> 28:16.600 very large, uh, are in thousands of dimensions. So, there are not many alternatives, uh, 28:16.600 --> 28:22.840 up to now to, to study those networks with, uh, sampling. So you have to use some, some, 28:23.080 --> 28:33.560 implementation like this. Okay, uh, just to wrap up, uh, this, uh, these algorithms for sampling 28:33.560 --> 28:41.400 and volume, uh, integration and implemented the C++. It's, uh, they're open source. So, 28:41.400 --> 28:50.040 germ scale is just the set, uh, in GitHub of, uh, all these, uh, packages, uh, it's the C++ 28:50.040 --> 28:57.560 package and the two interfaces. Uh, yes, I like these, uh, QR that the previous talks, uh, 28:57.560 --> 29:05.800 had, so I also had the QR. So, yes, this is supported by, uh, open source community 29:05.800 --> 29:13.800 and people that volunteer to, if, if they have, they find the interesting to, to work on this. 29:14.120 --> 29:22.840 Uh, so, this basic, basically how this is implemented and maintained. Okay, so, that's all. 29:22.840 --> 29:26.840 Thank you.