New Algorithms: Numerical Solution of Large Systems of Nonlinear Equations NSolve By Monodromy

Hello everyone. Thanks for joining this talk. I'm Arvind bhaskar, I'm a PhD candidate at university of Norde. I worked as a summer intern at Wolfram Research last summer and I worked on some new directions to numerical solution of large systems of nonlinear equations. I come at it from the perspective of engineering and applied mathematics. We try to solve some of these large systems in the context of kinematics of robotic systems, optimization problems for design of robotic mechanisms and so on. So this solver installed by Monodramy is a step in that direction and I'll be talking about the details of this solver. I created a Wolfram function repository by the same name with mentorship from Daniel Lichplau and Devendra Kaparia at Wolfram Research and I believe this is currently being integrated into NSOLVE in some form using some options. And if you had attended the morning talk by Devendra, he he mentioned about this in the upcoming version. So let me put that into some context. There are several methods to solve large nonlinear systems of equations. These could be root finding problem. They could also arise out of optimization problems because there is a duality there. Many global optimization problems can be posed as root finding problems of the gradient function to solve these equations. There is solve, there is n solve find root that has a whole plethora of methods. Algebraic elimination methods are available in solve and then eigen problems. Equivalent eigen problems are posed and solved in some unsolved methods. Find root is an option for local methods if you are looking at just a few instances of solutions to be found. We are interested in finding as many solutions as possible for these large systems. So that's why we started looking into some of these global methods that come under class of methods called numerical continuation methods. Some of these are available in NSOLVE as well some of these methods. But there are some new advancements in the class of numerical continuation methods that are more heuristic and incorporate a lot of probabilistic and randomness into the solution process and that gave birth to this new solver N solve by monodrami. It is iterative in nature. That's one of the distinguishing nature of these heuristic methods from other continuation methods. And we do not require construction of sophisticated start systems that are generally how classical numerical continuation methods work. And these can be applied to all sorts of systems that are complex analytic, not just restricted to polynomial systems. Let me first explain what I mean by monodramy because some some of you may not be familiar with it. This is an abstract Figure I mean, this comes from a live problem of finding cube roots of unity. And you need some abstraction. The circle at the bottom is what is a representation of the parametric space of this simple equation. Finding the roots of cube roots of unity. And the three colored balls you see here represent the three routes. Let's go to. Let's skip that slide. I welcome people to download this notebook and run it in their systems. I'll skip the slides where we have too much animation going on. And that was causing issues with the feed. Nice. Okay, so let's look at what kind of systems that we are studying. So we are looking at square systems and parametric equations. Because monodramy loops are always associated with some parameter space. You conduct these loops in the parameter space and study how the corresponding roots change as you do this and how these permutations evolve. And here I have given an algebraic description. F P is equal to 0, f is the system of equations square system where the dimension of F is same as the dimension of Z and P is a variable space associated with the system. To give an example, let's consider the equation of a circle shown here. Here X and Y represents the variable space and and P represents the parameter space and of this circle equation because you have the center of the circle and then the radius. So there are three parameters. Typically you conduct monotony loops in this parameter space, natural parameter space. But in some cases you might want to introduce artificial parameters and conduct monotonymy loops associated with the artificial parameter space. But when you do that in an artificial space, you lose some regularity assumptions and then that causes potentially issues for identifying a start root, we actually use find root function within Mathematica to get a start root. And that's how we do this. Depending on the parameter space, there are some associated risks. And typically we choose the constant parameter space because it is nice and well behaved. And then conduct monoramy loops using a circular construction as shown here. So these are complex. We conduct these monoramy loops in the complex parameter space because that that helps us with avoiding certain system singularities when the discriminant locus vanishes and so on. Other construction exists, but this is the one that we use. And for tracking these paths themselves, we use predictor corrector path tracking because these are high dimensional curve tracking. So starting from one root, you track a curve using predictor character tracking algorithms. There are sophisticated algorithms, but the key elements are two. One is you take a prediction based on the tangent information of the curve and then you correct it using a Newton Step I've given details of what, what techniques I use. Let me give a visualization here. Hopefully this, Hopefully this would not break. You can see I start from one root and then I do a curve tracking. I'm showing how the path tracking evolves and then you reach a new route. So imagine if you did not know where these routes lie, you just knew the one root starting route. By taking these random loops in the parameter space you pick, you can get to a new root and you don't have to stop at T is equal to one. You can produce, you can go further in the parameter space and then you get these nice crazy looking curves in the spirit of Halloween. And as you do track these curves, you actually get to many, many routes. I'll skip ahead again. I don't want that to crash my presentation again. The next thing we need is a root accumulation strategy. I've just shown you how you start from one root and then you can find many roots by randomly going around in the parameter space. Let's look at a concrete example here. I've shown a non algebraic equation. It has a pole at origin. So you can actually use this to solve even non holomorphic functions, some meromorphic functions with some poles in them. If people have complex analysis background. So let's say we want to find the roots of this system. We can already count here. 1, 2, 3, 4. There are these five roots within this domain shown here for this function. So let me. Run this. Anyway, so I'll add these illustrations. Be nice to show this illustration. Anyway, I'll skip ahead these slides. So we already. We also need some termination criteria because as these roots are new and newer and newer routes are completed, we need to figure out when to stop these things because we can keep newer and newer loops and finding more and more routes. So there are several criteria available in literature. Linear trace test that for polynomial systems that determines if we have found all. We could also use a simple criteria based on diminishing returns. Let's say you take further and further loops and then the number of routes that you are newer routes that you are getting diminishes. You can model those diminishing returns based off of a mark and recapture models. Lincoln Peterson estimates statistical models. These are the same models that are used in ecological studies where let's say you want to map out the number of species in a population, right? You take a sample of the species, tag them and release them back into the wild. And then after, let's say six months, you come back and capture some more and then see how many of them are tagged. And depending on the percentage of overlap you can project the total population in the. In the. In the neighborhood. Right. So same same thing we do with these root counts as well and figure out when to stop the algorithm. I'll get. And then we. Anyway, that's all the details. I've put all of that together. That's a lot of details there. Put that together into this NSL by Manorami solver. You can read about this if you go to Wolfram Function repository and look for this installed by Monodremy you will find this solver and some documentation as well. It has several options. The main option is this max roots criteria which is the novelty of this Using max criteria you can let's say for a large system, even with millions of routes, you want to just sample 100 roots. Let's say you can do that. It's very heuristic and our examples and the systems that we tested, fairly robust. But if somebody I know somebody express interest in knowing more about this in Devendra stock and if you use this and if it works, that's great. If it does not work, please get back to us and we'll see what can be done to address this. I'll show some examples where use this to effectively solve large systems of equations, not necessarily polynomial systems. This is an example of a 6R inverse kinematic problem. Let's say you are given this, you have this robot and you are given an end effector position and you want to figure out what are the inputs that would get you to this end effector position. This is a classical problem to a 11 dimensional system. And we show how NSOLVE by monotony can be used to solve for all the roots of this system in about four seconds. It has only 16 roots and we find them all. Nsolve can also find them, but in the current capability. But it returns a lot of us what are called spurious roots or repeated roots because it uses some other resultant based techniques that that result in. That result in spurious roots. And once you delete the duplicates you see that it corresponds to what ns albem monotony returns. And then you now I scale it up, take it to the limits, right? So this is an example of a kinematic problem that's well known in literature. Finding mechanisms that fit a set of points in plane. Let's say we are given nine points in plane and you want to design a mechanism that visits all these points. This is a geometric problem, but cannot be Solved easily, trivially. You can see the total algebraic degree of this system is crazy. This, you cannot solve this using some geometric constructions, obviously. And we show that for even though the algebraic degree is large, the actual number of roots is actually only 8,000. So classical numerical continuation methods do not scale well because you are constructing these large start systems and then a lot of the paths diverge and so on. But our heuristic continuation methods based in monadramy are able to find a large number of roots. In this case it. It is able to find 2000 roots in about 4300 seconds, however many minutes that is. And we see, we exhibit the results here. You can see a lot of these. These are all different solutions to the same problem found using N solve by monotony. So we found 2,000 roots for this system and it also extends naturally to transcendental systems. That's the one of the selling points of this, because you are not constructing any start system. You can actually do this to find roots or stationary points of multi dimensional systems. Let's see how this is. I hope this works. And this animation shows how the root generation generation evolves. So in this I'm showing one iteration of the monodomy loop. You can see all these roots are found starting from 1, 1 starting point. And then in the second second iteration, all of these roots proliferate in another direction and maps out all the stationary points. So we are actually exploiting the analyticity of the system. There are other ways to find these routes, like interval search, interval arithmetic methods. But you can see by exploiting analyticity, this scales better than interval arithmetic methods because interval arithmetic methods depend on they suffer from curves of dimensionality. Once the dimensions increases and you have to increase the search bound, it becomes really hard. And we show that even for a five dimensional transcendental system from literature, NSOL by Monarami can find as many as 45 roots in 26 seconds. But NSolve by itself, if you try to run it, it takes quite a bit of time and without returning any roots. So that's why integrating this capability in ensolve actually offers a lot more now. And I'm excited about hearing from the users and see how we can improve this further. And for my last part of the talk, I appreciate the patience. I'm sorry we had the technical difficulty. We show that this ensol by monotony and with some changes in heuristics, we actually solved an open research problem in literature that of finding all equilibrium configurations of a cable driven robot to give a brief Introduction to Cable Driven Robots Here you see an end effector platform that is operated by a lot of cables shown here and as you change the cable lengths the platform changes its configuration. But for the for a given set of cable lengths there can be multiple equilibrium configurations and it would be nice to know all of them so that we can design and plan these systems. And it's a 31 dimensional non algebraic system because of the cable sagging. The reason why it is non algebraic is because it has both algebraic and exponential parts. So I just shown a sub portion of the equations here. So alpha and beta if you see it occurs algebraically as well as in terms of hyperbolic sine and cosine functions. And we did better than state of the art solver using interval arithmetic and we increased the number of solutions significantly up to 50% in some of the test examples. I will add those animations to the I don't want to show them here because again I will crash the presentation. I'll add it to the handouts. I have also added a PDF file where you can learn more about this particular example. It's available as a preprint. The paper is currently under review. We have used these techniques to design and develop robotic systems for several other applications. Some of them I have shown here prosthetic application, there is repairs and so on. So I welcome you all to try this out and let me know your feedback and how it helped or did not help. I thank and acknowledge help from my managers advisors at Notre Dame as well as at Wolfram and I appreciate your patience and thanks for sticking.