Video 6-3: The simplex algorithm

0:04 Uh with that setup, I think uh I can tell you how the simplex algorithm works. Not in all full details in how it's implemented, but enough that you could I think work out the details if you wanted. And I realize some people have seen this before, maybe not in exactly that way. Um everyone explains it a little different. Uh I'm going to do a combination of of two explanations. first is um I'm trying to give you a geometric intuition because really LPS are not that hard. It's a very simple concept if you look at it uh from from this perspective. So here uh I'm going to take this first normal form where we have lesser equal constraints and potentially many of them u many more than variables and it still makes sense again here to make things simple um n is two so the number of variables is two I can still point it draw it in the plane each of these constraints um uh a lesser equal uh in in 2D is essentially one of those lines plus a direction uh which is not drawn but all of those kind of point in some way and to make the picture look sensible I guess you have to point them all in this way. So uh an individual constraint is basically telling you everything on this side is okay. If you take the intersection of all those half spaces then uh you get some some specific feasibility region. So here uh for notation I'm going to write hi for the half space of one constraint. So all the values or all vectors that satisfy this constraint. Okay. But really it's uh the area below a line or above a line in this case when n is two. Um and then the intersection of all these gives you this feasibility region the yellow part.

2:07 In this case again I'm I'm assuming that these all point towards the yellow bit uh which is not drawn in the picture. You can also see from the picture there's some constraints like in our Hessie James example that that don't even contribute to the feasibility region. They're redundant. Whereas um the actual contour of the feasibility region is formed by potentially many constraints. Okay. Now what about the cost function? In general, the cost function uh I drew it before as this dashed line and um it's in 2D it's some line with a certain slope but I can always visually at least think of rotating my picture so that this dashed line is a is a horizontal line. Okay. Now how to do the rotation with linear algebra math? I don't care for now. You can just uh change it to a different base. Um but suppose for now or maybe make it make it an assumption uh that the orange line is is horizontal which means if we move the objective value up and down, we move the orange line up and down which means the optimal value will be the lowest point in the feasibility region. And so you can physically imagine you put some marble in this yellow space and it'll just fall and roll down to the lowest point naturally. And here clearly it will end up here.

3:52 Okay. And again in 2D this is very simple but this actually works in general in in just the same way C itself. So uh brush up your linear algebra memory.

4:03 Um a vector in in full dimension defines a hyper plane namely all the um for example uh all the vectors that have uh scalar product zero with it and the same vector is also the normal vectors standing on on that plane and that's really what this green arrow is. So this is the um hyper plane for C and that's C itself. It's pointing downwards. So really C you can think of as a direction for improvement. That's how where you have to move and this works in in in any number of dimensions. So uh this is already showing a bit what the algorithm does, but uh imagine you have a marble that somehow sticks to the boundary. uh then it would still roll down step by step until it stops in the optimal point just following gravity or imagine a magnet pointing pointing downwards.

5:15 Okay. If you take um if you take enough half spaces um I guess this is hyper plane not half space. Half space is everything below a line. The hyper plane is the generalization of a line. And if you take um n such lines here too then you get a single point. Now what's the connection between halfspace and and hyper plane the halfspace has some boundary and that's a hyper plane here half space is some line below and then the line is the line boundary of the of the half space. So if you take n of these you always get a unique point. In this case, take any two uh uh constraints and you get the intersection point. Again, this also works in higher dimensions. And this intersection point is what we call a vertex. It could be outside of the feasibility region. It could be here or it could be one of the vertices of the feasibility region if you think of it as a polygon. Let me use the same same name. But really it's just the intersection of these uh hyperplanes. Um I guess that's that's a bit uh it's a bit wrong. These are the the hyper planes not the half spaces. They don't necessarily have just a single point at the intersection. Just from the picture uh from the physical intuition if you have some some point in the interior of this yellow yellow region you can always move it freely in all directions. So in particular you can move it in the direction of the improvement the green arrow. So you can always improve it further. So it's never going to be an optimal point. On the flip side that means you can always find an optimal point that is a vertex. The only other option would be that you have a whole edge that is optimal. Any point on that edge is optimal because you have a wall that's horizontal in that picture. Um but then still the end points of that wall would also be optimal. So in general you can say there's always an optimal vertex could be many or could be several. So the optimal objective value is given by the value of some vertex. So that's this x i star. And that gives you some Now we finally reduce this uh annoying continuous problem that I'm always a bit uh iffy about, not sure how to solve it into something discrete because there's only so many intersection points, namely as many ways as I have to uh as I can select um n hyperplanes from all the constraints. So these are all the vertices and in principle I could just go through all of them and if n is [Music] small that's actually poly time to brute force all the vertices is not how you want to solve this because in general n can be big I said millions in the beginning so certainly you don't want to do um m to the 1 million but uh If n is is really small like in our examples, it's one option. And that's what the simplex algorithm does in a more clever way. It doesn't brute force all the vertices, but rather it does it starts at a an arbitrary one and then always finds a neighbor that is better. And that means it follows some path along the contour and will stop as soon as there's no better neighbor. And then we're find we found an optimal vertex. Now what does a neighbor here mean? It means that it's defined by by an intersection of of hyperplanes where only one hyper plane is different. You throw out one and you add another in. Okay, that's the geometric view. And you can write this as an algorithm. um we're given uh the half spaces and then you can community say if if they're in a section is empty the problem is in feasible. Uh otherwise you find some feasible vertex. Uh how to do that is another round of simplex. Let's ignore that for now.

9:53 So suppose you're you're starting here with your X and then you look um if X can be improved or if the green direction is hitting a wall is going against some some constraint. If that's not the case then you can locally improve it. And so um if that's a funny bit um if all the uh neighbors that are feasible are actually strictly um strictly worse but you still can't locally improve then the problem is unbounded. That that case is a bit more tricky to see. Uh that needed another picture but um uh let's let's ignore that case for the intuition. The typical case is there is some way to move in the direction of the green line to uh a neighboring vertex and then you go there and repeat.

10:55 Okay. Now this is a nice way to visualize things and often times you get the right intuitive ideas. uh for example um we can't necessarily say how many points we visit if we follow such a contour and indeed this is in general a problem of simplex. You can find contrived examples where there is exponentially many vertices and you can also make simplex go to visit exponentially many vertices on its way to find the optimal solution. uh so often times the ge the geometry gives you the right uh way to think about the problem if you wanted to implement it it's not suitable uh for that you want to bring back linear algebra and so I want to briefly introduce that other side of things which is probably closer to what you've what you've seen if you've seen it uh the way um I'm going to present it is take the other normal form now with equality constraints that's a bit easier to work with because now really um you can work with these equalities and all that is uh is uh left of the other constraints is an in non- negativity constraint. That means we have uh now many more variables than constraints necessarily to make the problem interesting. If if the if that's not true, the problem is is essentially a either infeasible or has a single point or that that's easy to solve in all cases. So you should think now that there's many more variables than constraints. And now we talk about the U matrix A uh from linear algebra point of view. And for examp for for uh for starters we'll assume that it has full rank. So remember because M is now smaller the matrix A is is wide wider than tall uh because this is M and this is N. And so rank m means essentially that uh whenever you um or if it's non- degenerate and has full rank whenever you pick a collection of the of columns n time um m columns out of the n possible ones you get a set of linear linearly independent vectors. That's in general not true. But if you assume that the problem is such that uh you don't have redundant uh constraints or uh you know things are not there there's constraints that are either um multiples of each other or they intersect several in in a weird place.

13:40 If if you kind of ignore all these these corner cases then you can pick any collection of of columns and you get an a basis for the uh for the matrix which means just m independent linearly independent vectors. Now you can uh split up your variables. So again remember uh J is a collection of columns maybe these four um and then you can identify the corresponding variables that would put would be put on those columns as XJ. These are the basis variables, those that interact with these columns. And you can separate those from the others uh the non-base variables just by their index. Right? The the uh basis variables have an index in this capital J and the others don't. Now, uh why is that useful? um you can rewrite some of the things in splitting up the two parts uh to make progress you can split up the constraints. So again you take the matrix and now you uh take those four that gives your matrix AJ it's a square matrix now and all the others all the other columns together form this matrix AJ bar and you similarly do that for the constraint uh for the objective function coefficients then what we get is any feasible vector. So any any x that satisfies the equality constraints, you can just rewrite you split up uh that large sum into the columns that have index in capital J and the others. And so now this matrix AJ is square and has full rank.

15:57 which means it's invertible. So you can find an inverse matrix and you can multiply the equation from the left with this inverse matrix and then uh these cancel to identity and you can subtract this to the other side. So you get an expression for xj the basis variables as a function of the other variables. Okay. And that's of course uh it's kind of clear if you have uh a certain number of constraints but more variables you can fix some of the variables namely exactly uh m minus n many as that's really n minus m here right n minus m many and the other variables are then fully determined from the values you picked for the n minus m variables. Okay. Or why is that useful? Well, in particular, you can set all the other variables to zero. And then you get a relatively simple um simple uh vector x that solves the equality constraints that satisfies the equality constraints and these special kind of solutions. So the axis you get by taking some xj bar and setting it to zero and computing the others from this formula here uh that corresponds to the vertices in the geometric picture. Now they may or may not be a feasible solution because in general you don't know whether the base variables are satisfying the non- negativity constraints. They may be negative or not. But if if they are uh non- negative then they're a feasible basic solution.

17:47 So the wrap up to wrap up this slide uh if you're given a selection of columns any m columns out of the n variables you can compute the basic solution and you can check whether it's feasible just check if it's larger than zero. Now the last bit is uh to make this step in simplex check if we're optimal already and if not make an improving step and that um that is easy from where we are now. So we have this basic solution uh where we think of uh setting xj bar to zero. So this this whole term cancels. How can we locally modify this existing uh suppose it's feasible. How can we locally modify our solution to improve the objective function? Well, what we cannot do is easy to see. We uh we can't take um any of the uh xjs and and and modify them because we would violate the equality constraint. And uh we can't decrease the um non-base variables because that would violate the neg the non- negativity. They're already at zero. So the only thing is we can increase the non-base variables and then maybe see what happens to XJ if that makes it feasible or not. So we can rewrite the cost function and that will tell us what to do. So the that's the objective function. We can split it up the same way by the variable indices that occur in J and the others. Now we know XJ this is using for XJ what we uh what we know about the basic solution. So you plug this in in general. This is true not only for xj bar zero but in general any any solution that satisfies the equality constraints um has this form and now you can split up the two two terms. So you have this part which is this which is constant essentially there's no x in there and something um in front of xj bar and so this is the old cost vector and then this this new term here that combines uh sorry this shouldn't be here that combines uh these parts I guess the transpose is also Oh yeah, sorry. That should be just at the C. The matrix is already in the right shape. The the key point here is this is some vector CJ bar. Um I guess that's why the transpose. So this is transpose. Now uh you you get a a a new set of coefficients for an objective function. And uh the overall um the overall objective function okay I got to fix the I'm staying here on the corrected slide. Uh the overall objective function can be written as some constant plus some terms that depend only on the non-base variables. And these are the ones you're trying to um to modify. In particular, we try to slightly increase one of those and see what that does to the objective value.

21:16 Well, if the coefficient that some variable here has is is negative, then by increasing it, I'm only making okay, we minimize, then I'm making the uh objective function better. Whereas if the coefficient is positive, I cannot help the objective value by increasing this xj. So you get no local improvement exactly if this thing is is non- negative then any increase here would only make the cost higher and that's not better and indeed that is uh if you cannot locally improve is the same as having a global optimum. That's because everything here is convex. We have a convex function over a convex domain. So yeah, whenever you have a local optimum that's actually the global optimum. This is not true in general optimization. Uh but for l because of the linearity we we get this simple simple criterion. You're only looking for a local optimum. And so that's that's essentially uh the end of this simplex algorithm. That's what it does. It computes these CJs.

22:27 um it checks if if they're all greater or equal zero then the current solution is the optimum and you can stop otherwise there must be some entry that's negative and then you know which of these uh non-base variables now becomes a base variable and you have to kick out one of the base variables to move to a neighboring vertex and that's that means so you increase this variable and you do that um until sum of the base variables drops to zero and you can compute how far that uh that is. So you can essentially use more linear algebra to do that step in a smart way. But that's that's really getting into the implementation. So let me uh wrap up at this point. So what we've seen is um the conceptual ideas behind the simplex algorithm. Uh it's a relatively simple and and mostly combinatorial and discrete algorithm. it it uses linear algebra to do the computations. Um, but it's it's basically talking about these discrete objects as vertices on on a polyhedron. Polyhedron is the higher dimensional version of a poly of a polygon. It's fairly easy to implement once you dig through the the details what the computations have to do.

23:48 Implementing the implementations themselves is actually not hard. uh it's usually the the choice of open-source solvers um just because of that and it's also doing reasonably well anyways now what I've said before the worst case is actually exponential so simplex is not a poly time algorithm in the worst case uh the details depend on how you do this last step if you found some basis v some non-basis variable you want to increase and kick out another there's uh usually some choice and you have to pick a rule how to to make that choice and that depends um how long it takes but there's no such rule known that guarantees polomial time. On the other hand, it is known that if if you take any input and you slightly modify the numbers in it, you actually get an expected poly time algorithm with simplex. That's called smooth analysis. I'm just mentioning this here because it kind of fits the theme of of the lecture. uh but we're not going to dig into any of this in more detail. And last thing is there are alternative methods that are uh more complicated both to understand and to implement. Uh but they are achieving this worst case polomial time guarantee. Commercial solvers like Gurobi and so on they often implement those and these tend to perform best in large instances.

25:16 uh for the lecture all all you need to take from this is a a rough idea of what simplex does and the fact that you can solve LPS in polomial time in theory and in practice.