This work by Anthony Zhang is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
Copyright 2013-2017 Anthony Zhang.
Introduction to optimization.
Levent Tunçel
Section 081 (online)
Email: ltuncel@uwaterloo.caIn this course, we will be studying linear programming, integer programming, the Simplex algorithm, and the concept of duality, as well as various ways to efficiently approximate answers to these sorts of problems.
An abstract optimization problem (also known as a "P") is a problem in which we are trying to maximize or minimize a certain objective function - finding the inputs to the function, subject to certain constraints, such that the value of the function is maximised. In other words, given A \subseteq \mb{R}^n (the set of all feasible inputs, also known as the feasible region) and a function f: A \to \mb{R} (the objective function), find x \in A such that f(x) is the global maximum/minimum.
Optimization problems show up in almost every industry - hardware companies might want to simplify their circuit layouts, hotels might want to fill the most rooms possible at the highest price, and a factory might want to generate the most profit given constraints on materials and labour. Also, you can use it to play Pandemic optimally.
There are three types of abstract optimization problems:
Abstract optimization problems are usually given as prose descriptions of the problems, which we then turn into abstract mathematical models/programs/problems. Then, we can take advantage of the many powerful computer solvers available today to automatically find optimal solutions.
Suppose a company produces either plastic carrots (which sell for $5 each and take 4 minutes to make) or rubber potatoes (which sell for $8 each and take 8 minutes to make). The company has 500 minutes available in the rest of the day, and the machine that makes either of these products can only be used up to 100 times - how many of each product should be made to maximize profit?
First, we define our decision variables - the actual values we are trying to find. Let c be the number of carrots, and p be the number of potatoes.
Profit can now be defined as a function of those decision variables: f(c, p) = 5c + 8p.
Next, we can define constraints to specify A. There are only 500 minutes left, and the time spent depends on the number of each product made, so 4c + 8p \le 500. There are only 100 products left that we can make, and this is simply the number of carrots and potatoes, so c + p \le 100. Additionally, c and p must be non-negative, since they represent the number of each product.
We now have a mathematical model of the problem: "maximize 5c + 3p subject to 4c + 3p \le 500 and c + p \le 100".
A feasible solution to an abstract optimization problem is any x \in A - any variable assignment such that the constraints of the problem are satisfied, like c = 50, p = 10. An optimal solution is one that actually solves the problem - a feasible solution that results in the global maximum/minimum of the objective function. Likewise, feasible solutions to prose descriptions of the problem are assignments to the unknowns, like "make 50 carrots and 10 potatoes".
To show that a mathematical model of a problem is correct (properly represents the prose description of the problem), we usually define a bijection between feasible solutions for the mathematical model and feasible solutions for the prose description of the problem, such that the bijection preserves cost - each x \in A for the mathematical model must have f(x) equal to the the profit in the prose description for its corresponding feasible solution.
For example, show that the model from the above example is correct:
4c + 3p \le 500 \land c + p \le 100 can be trivially bijected to "carrots take 4 minutes to make and potatoes take 8 minutes, and the machine can only be used to make 100 products".
Clearly, f(c, p) = 5c + 8p has the same value as the profit since "carrots sell for $5 each and potatoes sell for $8 each".
In this course we will mostly focus on affine optimization problems, since solving general optimization problems is much harder.
An affine function is one that can be written as f(\vec{x}) = \vec{a} \cdot \vec{x} + \beta where \vec{a} \in \mb{R}^n and \beta \in \mb{R}. If \beta = 0, the function is also linear.
A linear program is an abstract optimization problem of the form \min \set{f(x) : x \in \mb{R}^n \wedge \left(\forall 1 \le i \le m, g_i(x) \le b_i\right)} where f is affine, g_1, \ldots, g_m are all linear, and m is finite. We might also replace \min with \max depending on the problem.
Usually we write this as: "minimize f(x) subject to g_1 \le b_1, \ldots, g_m \le b_m", or "maximize f(x) s.t. g_1 \le b_1, \ldots, g_m \le b_m". It's also common to use \vec{x} \le \mb{0} as shorthand for x_1 \ge 0, \ldots, x_n \ge 0.
Suppose an oil company in the next 4 months needs to supply 5000L, 8000L, 9000L, and 6000L, respectively, and oil costs them 0.75/L, 0.72/L, 0.92/L, and 0.90/L, respectively. The company has an oil tank that can store up to 4000L, and starts off with 2000L. Oil is bought at the beginning of the month, is used throughout the month, and is all dumped into the tank at the end of the month. What quantity of oil should be bought in each month to minimize cost?
Clearly, the inputs are the quantity of oil bought in each month, b_1, b_2, b_3, b_4. For convenience, we will also define the tank value at the beginning of each month, t_1, t_2, t_3, t_4. Clearly, 0 \le t_1 \le 4000, 0 \le t_2 \le 4000, 0 \le t_3 \le 4000, 0 \le t_4 \le 4000, since the capacity of the tank is 4000L.
The objective function we are minimizing is f(b_1, b_2, b_3, b_4) = 0.75b_1 + 0.72b_2 + 0.92b_3 + 0.90b_4.
Clearly, the tank level in month i is t_1 = 2000, t_2 = t_1 + b_1 - 5000, t_3 = t_2 + b_2 - 8000, t_4 = t_3 + b_3 - 9000, since the tank level is just the leftover oil after supplying the last month's oil.
Clearly, in order to supply enough oil, we must have t_1 + b_1 \ge 5000, t_2 + b_2 \ge 8000, t_3 + b_3 \ge 9000, t_4 + b_4 \ge 6000. However, the first three constraints are redundant with the previously defined constraints t_2 = t_1 + b_1 - 5000, t_3 = t_2 + b_2 - 8000, t_4 = t_3 + b_3 - 9000 and 0 \le t_1, 0 \le t_2, 0 \le t_3, 0 \le t_4, so we only need the t_4 + b_4 \ge 6000 constraint.
These encode all of the problem constraints, so we have "minimize 0.75b_1 + 0.72b_2 + 0.92b_3 + 0.90b_4 subject to 0 \le t_1 \le 4000, 0 \le t_2 \le 4000, 0 \le t_3 \le 4000, 0 \le t_4 \le 4000 and t_1 = 2000, t_2 = t_1 + b_1 - 5000, t_3 = t_2 + b_2 - 8000, t_4 = t_3 + b_3 - 9000 and t_4 \ge 6000".
If we use a linear programming solver, we get b_1 = 3000, b_2 = 12000, b_3 = 5000, b_4 = 6000, which gives us the minimum cost 20890.
Same thing, but minimise the maximum amount of oil purchased in any single month, and then show that the model is correct assuming that the above example's model was correct:
From the previous example, we have the same constraints 0 \le t_1 \le 4000, 0 \le t_2 \le 4000, 0 \le t_3 \le 4000, 0 \le t_4 \le 4000 and t_1 = 2000, t_2 = t_1 + b_1 - 5000, t_3 = t_2 + b_2 - 8000, t_4 = t_3 + b_3 - 9000 and t_4 \ge 6000.
However, the objective function is different - let M be the maximum amount of oil bought in a single month, and we can minimize this value. We can define M to be the maximum amount of oil bought in a single month by adding the constraints b_1 \le M, b_2 \le M, b_3 \le M, b_4 \le M.
These encode all of the problem constraints, so we have "minimize M subject to 0 \le t_1 \le 4000, 0 \le t_2 \le 4000, 0 \le t_3 \le 4000, 0 \le t_4 \le 4000 and t_1 = 2000, t_2 = t_1 + b_1 - 5000, t_3 = t_2 + b_2 - 8000, t_4 = t_3 + b_3 - 9000 and t_4 \ge 6000 and b_1 \le M, b_2 \le M, b_3 \le M, b_4 \le M".
Now to show correctness. Suppose that M is an optimal solution that has values b_1, b_2, b_3, b_4. Since M is a feasible solution, M \ge \max\set{b_1, b_2, b_3, b_4}, and since M is an optimal solution, it is the lowest such value that satisfies that, so M = \max\set{b_1, b_2, b_3, b_4}. Therefore, M is also the minimum maximum amount of oil purchased in any single month.
Integer programs are just linear programs plus a integrality constraint - a constraint that all of the variables must be integers. When we formulate integer programs, we just add the constraint x_1, \ldots, x_n \text{ are integral}.
An integer program is mixed if there are still non-integral variables, and pure if it only has integral variables. Integer programs are harder to solve than LPs.
The size of a linear/integer program generally refers to either the number of variables, of the number of constraints. The running time of an IP/LP solver algorithm is the number of steps it takes, as a function of the size of the IP/LP.
Efficient algorithms are polynomially. All LPs can be solved in polynomial time, though with a high exponent. Solving integer programs, however, are in NP, and so are highly unlikely to be solvable in polynomial time.
The knapsack problem is a special case of integer programs. Given a set of n items, with values v_1, \ldots, v_n and weights w_1, \ldots, w_n, how do we maximize the total value of the items we can carry in a knapsack given that the total weight cannot exceed m? From CS341, we know that if the number of items can be fractional, we can just use a greedy algorithm.
As an integer program, we have: maximize v_1 x_1 + \ldots = v_n x_n subject to w_1 x_1 + \ldots + w_n x_n \le m, x_1, \ldots, x_n \text{ are integral}.
Write the IP for a knapsack problem where we have n = 4 and either x_1 + x_2 is at least 4, or x_3 + x_4 is at least 4, and we can only send x_3 if we send at least one x_4:
We're trying to maximize v_1 x_1 + v_2 x_2 + v_3 x_3 + v_4 x_4.
Clearly, we can express x_1 + x_2 \ge 4 AND x_3 = x_4 \ge 4, but how do we express x_1 + x_2 \ge 4 OR x_3 + x_4 \ge 4?
One way to do this is by adding a binary variable y to our IP, constrained such that 0 \le y \le 1 and y \text{ is integral}.
We can now express x_1 + x_2 \ge 4 \lor x_3 + x_4 \ge 4 as x_1 + x_2 \ge 4y, x_3 + x_4 \ge 4(1 - y). Whenever y is 0, 4y is 0, so x_1 + x_2 will always be true (since x_1 and x_2 are non-negative), while 4(1 - y) is 4, so the second constraint becomes x_3 + x_4 \ge 4. Whenever y is 1, 4y is 4, so the first constraint becomes x_1 + x_2 \ge 4, while the second constraint is always true (since x_3 and x_4 are non-negative).
Now, how do we express sending x_3 only if we send at least one x_4? We want to express that x_4 = 0 \implies x_3 = 0.
One way to do this is to multiply x_4 by a large enough value so that if it's 1 or more, the resulting value will definitely be greater than x_3. For example, x_3 must be no greater than \floor{\frac m {w_3}}, so we can add the constraint x_3 \le \floor{\frac m {w_3}} x_4 to do this. When x_4 = 0, \floor{\frac m {w_3}} x_4 is also 0, so x_3 must be 0 as well. When x_4 \ge 1, \floor{\frac m {w_3}} x_4 \ge \floor{\frac m {w_3}}, so x_3 \le \floor{\frac m {w_3}} x_4 is always true.
Therefore, we can write this as "maximize v_1 x_1 + v_2 x_2 + v_3 x_3 + v_4 x_4 subject to w_1 x_1 + w_2 x_2 + w_3 x_3 + w_4 x_4 \le m, x_1 + x_2 \ge 4y, x_3 + x_4 \ge 4(1 - y), 0 \le y \le 1, y \text{ is integral}, x_3 \le \floor{\frac m {w_3}} x_4, x_1, x_2, x_3, x_4 \ge 0, x_1, x_2, x_3, x_4 \text{ is integral}".
The above example demonstrates a common trick of using a binary variable to express logical constraints involving disjunction. The binary variable lets us make constraints that apply when it has one value, and always be true for other values.
Scheduling is one of the most common optimization problems in industry. Suppose we have the coffee shop that is open on weekdays, hires workers that work for 4 consecutive weekdays per week (can wrap around weeks, and ignoring weekends), and requires 3, 5, 9, 2, and 7 workers to meet demand on Monday through Friday, respectively. How do we hire the smallest number of workers that meets this demand, and how do we schedule them to meet this demand?
First, we could define 5 variables, x_1, \ldots, x_5, for the number of workers who should start working on each weekday. Then, the objective becomes to minimize x_1 + \ldots + x_5.
What are the feasible solutions? Well, the constraints are that the number of workers starting on each day must be non-negative, and that the demand is met each day. To ensure we have enough demand for each day, we can define five constraints for each weekday. For example, for Monday demand to be met, we need the number of workers who start on Monday, Friday, Thursday, or Wednesday to be at least 3, since Monday demand is 3 and on Monday, we have workers from up to 4 days ago: x_1 + x_5 + x_4 + x_3 \ge 3.
So, the problem becomes: "minimize x_1 + \ldots + x_5 subject to x_1 + x_5 + x_4 + x_3 \ge 3, x_2 + x_1 + x_5 + x_4 \ge 5, x_3 + x_2 + x_1 + x_5 \ge 9, x_4 + x_3 + x_2 + x_1 \ge 2, x_5 + x_4 + x_3 + x_2 \ge 7, and x_1, \ldots, x_5 \ge 0".
Suppose we wanted to constrain a variable x to be one of k_1, \ldots, k_n. One way we could do this is to have binary variables y_1, \ldots, y_n, where y_i represents x = k_i - adding a constraint of the form x = k_1 y_1 + \ldots + k_n y_n. Then, we can ensure that exactly one of them is true by adding the constraint y_1 + \ldots y_n = 1.
Another common type of linear programming problem in the real world is graph optimization. For example, how do we find the shortest path on a map? How do we optimize a microchip layout?
For the shortest path on a map example, we can associate street intersections with vertices, and streets with edges, weighted by the distance (a non-negative real number) along its length between intersections.
An s, t-walk in a graph is a sequence of edges s v_1, v_1 v_2, \ldots, v_k t. An s, t-path is an s, t-walk such that non-consecutive edges in the sequence don't share any vertices (we never revisit the same vertex twice).
The length of an s, t-path is the sum of the weights of the edges in the path. So, given a graph G = \tup{V, E}, a weight function w: E \to \mb{R}, and two vertices s, t, we want to find the minimum length s, t-path. How do we write this as an IP?
Suppose we have tasks t = t_1, \ldots, t_n and machines m = m_1, \ldots, m_n, and c_{m_i, t_j} represents the cost of having machine m_i perform task t_j. How do we match machines to tasks to minimize total cost?
Create a graph with one vertex for each job and edges between them weighted by the corresponding cost.
We want to find the perfect (every vertex in the graph is incident to an edge in the matching) matching (subset of edges such that none share vertices) in the bipartite graph with the minimum total weight.
Let \delta(v) be a function that returns the set of all edges incident to a given vertex. Then a set of edges M is a perfect matching if and only if \forall v \in V, \magn{M \cap \delta(v)} = 1.
For our IP, we can have a binary variable x_{\text{machine}, \text{task}} for every edge \tup{\text{machine}, \text{task}}. Now we can simply minimize \sum_{m_i \in m, t_j \in t} c_{m_i, t_j} x_{m_i, t_j}.
To constrain solutions to perfect matchings only, we can add the constraints \sum_{m_i \in m} \sum_{t_j \in \delta(m_i)} x_{m_i, t_j} = 1 and \sum_{t_j \in t} \sum_{m_i \in \delta(t_j)} x_{m_i, t_j} = 1 (for each vertex, choose exactly one incident edge).
The shortest path problem: given a graph G = \tup{V, E}, what is the s, t-path that has the minimum total length? We want to formulate this as an IP.
Recall from MATH239 that a cut in a graph is a partitioning of V into two disjoint subsets, often written as just one of the sets. In other words, a cut is a pair \tup{S, T} such that S \subseteq V, T \subseteq V, S \cap T = \emptyset, S \cup T = V.
An s, t-cut is a set of edges rather than a pair of sets of vertices. Let \delta(S) = \set{\set{u, v} \in E : u \in S, v \notin S} for S \subseteq V. Then given s \in S and t \notin S, \delta(S) is an s, t-cut. To enumerate all s, t-cuts, we can simply find all subsets of V that contain s but not t, then apply \delta to those sets of vertices.
Clearly, if we remove the edges contained in the s, t-cut from the graph, s and t would no longer be connected. Therefore, every s, t-path must contain at least one edge from every s, t-cut.
Also, if S \subseteq E contains at least one edge from every s, t-cut, then S also is a superset of an s, t-path. This is because if S had an edge from every s, t-cut and didn't have an s, t-path, the set of all vertices R that can be reached from s is an s, t-cut (since it includes s but not t, as no s, t-path exists by assumption), so S must include an edge from \delta(R). However, S can't contain any edges from \delta(R), because if it did contain an edge \set{u, v} then both u and v should have been in R (since they're reachable from s). Therefore, S must also contain an s, t-path.
Define binary variables x_1, \ldots, x_m for each edge e_1, \ldots, e_m being or not being in the shortest path, respectively. Let w_1, \ldots, w_n be the weights of each edge e_1, \ldots, e_m, respectively.
Clearly, we want to minimize x_1 w_1 + \ldots + x_m w_m (i.e., the length of the path). What are the feasible solutions (i.e., supersets of s, t-paths)?
To be an s, t-path, the edges of the path must include at least one edge from every s, t-cut. This is a lot easier to write as constraints: we have one constraint for each s, t-cut \delta(S), and each constraint ensures that at least one of the edge variables T \subseteq \set{x_1, \ldots, x_m} corresponding to \delta(S) is true.
For example, given the graph G = \tup{\set{a, b, s, t}, \set{ab, as, at, bs, bt}} with weights w_1, \ldots, w_5, we might write the shortest path problem as "minimize x_1 w_1 + \ldots + x_5 w_5 subject to x_2 + x_4 \ge 1 (for the s, t-cut \set{sa, sb}), x_1 + x_2 + x_3 \ge 1 (for the s, t-cut \set{ab, as, at}), x_1 + x_2 + x_5 \ge 1 (for the s, t-cut \set{ab, as, bt}), x_3 + x_5 (for the s, t-cut \set{at, bt}), and x_1, \ldots, x_m are non-negative integers".
Note that we don't need an upper bound on x_1, \ldots, x_m. This is because if it was greater than zero, then it wouldn't be an optimal solution.
Nonlinear programs are problems of the form "minimize f(x) subject to g_1(x) \le 0, \ldots, g_m(x) \le 0" where f, g_1, \ldots, g_m: \mb{R}^n \to \mb{R}. This is a superset of linear programs and integer programs. Basically, we can do arbitrary constraints, math permitting.
For example, we can minimize Euclidean distance by using an objective function (x_1 - p_1)^2 + (x_2 - p_2)^2 + (x_3 - p_3)^2. Another example is constraining variables to integers: simply add the constraint \sin(\pi x) = 0 to ensure that x is an integer, or x(1 - x) = 0 to ensure that x is either 0 or 1.
It's even possible to write Fermat's last theorem as a single NLP: "minimize (a^n + b^n - c^n)^2 subject to \sin(\pi a) = 0, \sin(\pi b) = 0, \sin(\pi c) = 0, \sin(\pi n) = 0, a, b, c \ge 1, n \ge 3" - the objective function solution is 0 if and only if a^n + b^n = c^n, and the theorem is true if and only if the optimal value is greater than 0.
What does it mean to solve an optimization problem? What is a solution to an optimization problem?
Consider the LP "maximize 2x_1 + 3x_2 subject to x_1 + x_2 \le 1, x_1, x_2 \ge 0". Clearly, in this case the optimal solution is simply x_1 = 1, x_2 = 0. However, in general, it is more difficult to say what exactly can be considered a solution.
A feasible solution is an assignment of values to each of the variables in the LP such that all of the constraints are satisfied - regardless of the objective function value. A feasible optimization problem is one that has at least one feasible solution, and an infeasible optimization problem is one that has none. An optimal solution is a feasible solution where the objective function is maximized (for a maximization problem), or minimized (for a minimization problem).
However, an optimization problem having a feasible solution doesn't mean that it has an optimal solution - the maximum/minimum objective function value could be unbounded, like for "maximize x subject to x \ge 0". If for any feasible solution there exists another with a higher/lower objective function value (depending on whether it's a maximization/minimization problem), the optimization problem is unbounded, and otherwise, it is bounded.
Consider the optimization problem "maximize x subject to x < 1". Clearly, this is feasible (x = 0 is a solution), and bounded (x \le 1 because x < 1), yet it has no optimal solution - the objective function can always get closer to 1 for any given value of x! The same occurs for the optimization problem "minimize \frac 1 x subject to x \ge 1". Note that neither of these are LPs, because one uses a strict inequality, and the other a division.
LPs are a bit better behaved than general optimization problems. According to the Fundamental Theorem of Linear Programming, an LP either has at least one optimal solution, is infeasible, or is unbounded.
An LP solver should therefore output one of the following: an optimal solution to the input LP (plus proof that it's optimal), a proof that the LP is infeasible, or a proof that the LP is unbounded.
Consider the linear system A\vec x = \vec b, \vec x \ge \vec 0. If we can find a \vec y such that {\vec y}^T A \ge {\vec 0}^T and {\vec y}^T \vec b < {\vec 0}^T, then the system can't have any solutions. In other words, we can add/subtract/scale certain equality constraints in the linear system together to get a new constraint such that the coefficients of the constraint are all non-negative, yet the right hand side is negative. Here, \vec y simply represents the fact that we're adding/subtracting/scaling rows in A and elements in \vec b to form a new constraint.
If we have such a constraint (left side coefficients non-negative, right side negative), then the left side of that constraint must be a non-negative value (since the elements of \vec x are also non-negative) and the right side is negative (by assumption), so this constraint is not satisfiable - the LP must be infeasible.
According to Farkas' Lemma, if there's no feasible solutions, then a \vec y that satisfies those conditions must exist. Therefore, a value of \vec y such that {\vec y}^T A \ge {\vec 0}^T and {\vec y}^T \vec b < 0 is a proof of infeasibility. For an LP in SEF, a value \vec y such that {\vec y}^T A = {\vec 0}^T and {\vec y}^T \vec b \ne 0 works just as well for proof purposes - a linear combination of constraints that results in the left side being zero and the right side nonzero.
For example, consider an LP with constraints 3x_1 - 2x_2 = 6, 2x_1 - x_2 = 2. If we choose y = \begin{bmatrix} -1 \\ 2 \end{bmatrix}, then we get {\vec y}^T A = \begin{bmatrix} 1 & 0 \end{bmatrix} and {\vec y}^T \vec b = -4. This means that the two constraints, with the first one subtracted from the second one doubled, imply the constraint x_1 = -2, which is always infeasible because \vec x \ge \vec 0.
To prove that an optimal solution \vec x for an LP is actually optimal, we simply need to prove that for any feasible solution \vec x', the objective function value is less or equal to the objective function value for \vec x - basically, proving that the upper bound of the objective function occurs at \vec x. This is a proof of optimality. We will look at how to construct these proofs systematically later on, using the strong duality theorem.
Consider the LP "maximize 2x_1 + 3x_2 subject to x_1 + x_2 \le 1, x_1, x_2 \ge 0". Clearly, any feasible solution x_1, x_2 must be
To prove that an LP is unbounded, we construct a family of feasible solutions \vec x(t) for all t \ge 0, then show that as t goes to infinity, so does the value of the objective function. In other words, we need to prove that for any arbitrarily high objective function value (assuming a maximization problem), we can construct a variable assignment that is both a feasible solution, and results in that objective function value - a proof of unboundedness.
This family of feasible solutions is always a line in the linear space. Therefore, a family of feasible solutions \vec x(t) = \vec x_0 + t \vec r such that for any z, there exists a t such that \vec x(t) = z is a proof of unboundedness. In matrix form, a proof of unboundedness is simply a value of \vec x_0, \vec r such that \vec x_0 \ge \vec 0, \vec r \ge \vec 0, A \vec x_0 = \vec b, A \vec r = \vec 0, {\vec c}^T \vec r > 0. Also, if the LP is unbounded, those values of \vec x_0, \vec r must exist.
Consider the LP "maximize x subject to x \ge 0". Clearly, for any objective function value z \ge 0, we can construct an assignment x = z, which must be a feasible solution since z \ge 0. Additionally, the objective function value for the feasible solution x must be z. Therefore, the LP is unbounded, and a proof for that would be \vec x_0 = \begin{bmatrix} 0 \end{bmatrix}, \vec r = \begin{bmatrix} 1 \end{bmatrix}.
Basically, whichever outcome an LP has, there must exist a proof that that outcome is the case.
A free variable is a variable that is not constrained.
An LP is in standard equality form (SEF) if and only if all of the following are true:
In other words, an LP in SEF is a problem of the form "maximize \vec c \cdot \vec x subject to A \vec x = \vec b, \vec x \ge \vec 0", where A is an n by m coefficients matrix and \vec x, \vec b, \vec c \in \mb{R}^m.
Every LP has an "equivalent" LP that is in SEF. By equivalent, we mean the original LP is unbounded if and only if the equivalent LP is unbounded, the original LP is infeasible if and only if the equivalent LP is infeasible, and an optimal solution of one can easily be used to compute an optimal solution of another.
To convert an LP into its equivalent in SEF:
Once this is done, we have an equivalent LP in SEF. Also, we can easily convert optimal solutions for the LP in SEF into optimal solutions for the original LP, simply by taking the values of the corresponding variables and ignoring the newly introduced variables.
We can actually solve LPs by hand, by building an intuition for what would maximise the objective function, and trying to satisfy the constraints with that in mind.
For example, suppose we have the LP "maximise 4x_1 + 3x_2 + 7 subject to 3x_1 + 2x_2 + x_3 = 2 and x_1 + x_2 + x_4 = 1 and x_1, x_2, x_3, x_4 \ge 0". Clearly, x_1 = 0, x_2 = 0, x_3 = 2, x_4 = 1 is a feasible solution, because if we set x_1 = x_2 = 0, the constraints become x_2 = 2, x_4 = 1. This has objective value 7. How do we find a feasible solution with a higher objective function value?
To get a better objective function value, we need to either increase x_1 or x_2, the two variables in the objective function. To keep things simple, we'll just arbitrarily choose to modify only one of the variables for now - maximize x_1 while keeping x_2 = 0.
Now the LP becomes "maximise 4x_1 + 7 subject to 3x_1 + x_3 = 2 and x_1 + x_4 = 1 and x_1, x_3, x_4 \ge 0". To maximize this, we just need to maximize x_1. Clearly, the LP is feasible whenever x_3 = 2 - 3x_1 and x_4 = 1 - x_1 - we have a feasible solution for every possible value of x_1 where all the non-negativity constraints are satisfied.
Since x_3 \ge 0 and x_3 = 2 - 3x_1, 2 - 3x_1 \ge 0, so x_1 \le \frac 2 3. Since x_4 \ge 0 and x_4 = 1 - x_1, 1 - x_1 \ge 0, so x_1 \le 1. Therefore, since we also know that x_1 \ge 0, 0 \le x_1 \le \frac 2 3. Since we're trying to find the largest value of x_1 that still satisfies the constraints, x_1 = \frac 2 3. This gives us an objective function value of \frac 8 3 + 7 with the solution x_1 = \frac 2 3, x_2 = 0, x_3 = 0, x_4 = \frac 1 3, which is better than our original solution.
Turns out, however, this solution is still not optimal. Additionally, we can't increase x_2 without decreasing x_1 at this point, so we can't just repeat the trick with x_2 this time instead of x_1. However, it would be possible to use the "maximize one variable while keeping others in the objective function the same" trick again if we rewrote the LP in what's known as canonical form. We will look at how to do this later.
The basic idea behind the simplex algorithm:
The main questions are: What's canonical form and how do we rewrite LPs in that form? Why does it ensure that we can find a better feasible solution?
The column sub-matrix notation allows us to select a subset of columns from a matrix into a new matrix: if A = \begin{bmatrix} a_{1, 1} & \ldots & a_{1, m} \\ \vdots & \vdots & \vdots \\ a_{n, 1} & \ldots & a_{n, m} \end{bmatrix} and B = \set{b_1, \ldots, b_k}, 1 \le b_i \le m, then A_B = \begin{bmatrix} a_{1, b_1} & \ldots & a_{1, b_k} \\ \vdots & \vdots & \vdots \\ a_{n, b_1} & \ldots & a_{n, b_k} \end{bmatrix}. As a short form, we can also do A_j where 1 \le j \le m, which means "column j of the matrix A", a vector of size n.
For a vector \vec x = \begin{bmatrix} x_1 \\ \vdots \\ x_n \end{bmatrix}, we also often use the subscript to denote the row sub-matrix \vec x_B = \begin{bmatrix} x_{b_1} \\ \vdots \\ x_{b_k} \end{bmatrix}. However, this isn't that common and it's generally preferred to write it out in full to avoid confusion.
A basis of a matrix is a subset of column indices B \subseteq \set{1, \ldots, m} such that A_B is a square matrix and A_B is non-singular. In other words, B is a set of size n, and the columns of A_B are linearly independent (so they can't be written as linear combinations of each other). Not all matrices will have a basis; for example, a 2 by 2 zero matrix (the only value of B that could result in a square A_B is \set{1, 2}, but A_B is singular).
We often denote the inverse of the basis as \overline B = \set{1, \ldots, m} - B - the column indices for columns that are not in the basis B.
From MATH136 we learned that the number of linearly independent columns in a matrix is the same as the number of linearly independent rows. This means that the number of independent columns is between 0 and the number of independent rows. Therefore, if we have a square matrix A_B, it can only possibly be a basis if all of the rows of A are independent; if there were less than n independent rows, there would also be less than n independent columns, so A_B must be singular. Therefore, a matrix has a basis if and only if all of its rows are independent, and that basis is a set of n independent columns in the matrix (which is a maximal set of independent columns in the matrix).
Suppose we have an LP in standard equality form: "maximize {\vec c}^T \vec x subject to A \vec x = \vec b". If B is a basis of A, then a variable x_i is a basic variable for B if i \in B, or a non-basic variable for B if i \notin B (i.e., i \in \overline B). Essentially, the coefficients of each non-basic variable can be written as a linear combination of the coefficients of the basic variables.
A basic solution for an LP and a basis B is a solution \vec x for the LP such that A \vec x = \vec b and all non-basic variables are set to 0. For the basic solution for a basis B = \set{b_1, \ldots, b_k}, the LP becomes "maximize {\vec c}^T \vec x subject to A_B \begin{bmatrix} x_{b_1} \\ \vdots \\ x_{b_k} \end{bmatrix} = \vec b".
Since B is a basis, A_B must be invertible, so A_B^{-1} A_B \begin{bmatrix} x_{b_1} \\ \vdots \\ x_{b_k} \end{bmatrix} = A_B^{-1} \vec b, or \vec x_B = A_B^{-1} \vec b. Interestingly, this means that a basis B has exactly one basic solution, and that solution is a vector \vec x where {\vec x}_B = A_B^{-1} \vec b and {\vec x}_{\overline B} = \vec 0.
Note that a basic solution doesn't necessarily satisfy the non-negativity constraints \vec x \ge \vec 0. A basic solution \vec x is feasible if and only if it does satisfy this constraint \vec x \ge \vec 0, because by definition it already satisfies A \vec x = \vec b. Solutions can be any combination of feasible/not-feasible and basic/not-basic.
A basic solution \vec x for an LP (without a specified basis) is a feasible solution that is a basic solution for the LP for one or more bases for that LP. So for some set of n linearly independent columns \vec x has corresponding values that result in the constraints being satisfied, and for all other columns \vec x has the element 0.
Suppose again that we have an LP in standard equality form: "maximize {\vec c}^T \vec x subject to A \vec x = \vec b". If the rows of A are dependent (at least one row can be written as a linear combination of others), then either the LP is infeasible, or A has a redundant constraint, which we can remove. Here's why: since there are dependent rows, one of the rows could be written as a linear combination of others, so we could cancel out its coefficients by adding/subtracting/scaling constraints to get the new constraint 0x_1 + \ldots + 0x_m = k. If k \ne 0, then the constraint is clearly not satisfiable, so the LP is infeasible, and if k = 0, then that constraint is redundant because it just says 0 = 0.
Suppose again that we have an LP in standard equality form: "maximize {\vec c}^T \vec x subject to A \vec x = \vec b". Suppose B is a basis for A.
The LP is in canonical form if and only if A_B = I (the column sub-matrix is the identity matrix) and c_{b_1} = \ldots = c_{b_k} = 0 (the objective function only contains non-basic variables).
Since A_B = I for an LP in canonical form with respect to a basis B, A_B^{-1} = I^{-1} = I, so the basic solution for B is just \vec x such that {\vec x}_B = \vec b and {\vec x}_{\overline B} = \vec 0.
Interestingly, it's always possible to write an LP into canonical form given a basis, such that the new LP has the same feasible region and feasible solutions have the same objective function value as the original LP. How do we do this?
Suppose we have the LP in SEF "maximize \begin{bmatrix} 0 & 0 & 2 & 4 \end{bmatrix} \vec x subject to \begin{bmatrix} 1 & 0 & 1 & -1 \\ 0 & 1 & 1 & 2 \end{bmatrix} \vec x = \begin{bmatrix} 1 \\ 2 \end{bmatrix} and \vec x \ge \vec 0". One basis for this LP is B = \set{2, 3} - it's a square matrix, and its columns are linearly independent.
To convert this LP into canonical form given B, we start by making A_B into the identity matrix. To do this, we subtract the first equality constraint from the second, and then swap them to get the equality constraints \begin{bmatrix} -1 & 1 & 0 & 3 \\ 1 & 0 & 1 & -1 \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}.
In general, to do this we left multiply both sides of A \vec x = \vec b by A_B^{-1} to get A_B^{-1} A \vec x = A_B^{-1} \vec b. Since A_B^{-1} A is the matrix formed by multiplying A_B^{-1} with the columns of A individually, the columns in B will end up as columns of the identity matrix, so (A_B^{-1} A)_B is indeed the identity matrix. Because A_B is non-singular, A_B^{-1} A \vec x = A_B^{-1} \vec b has the exact same solutions for \vec x as A \vec x = \vec b, too.
Now, we need to make \begin{bmatrix} c_{b_1} \\ \vdots \\ c_{b_k} \end{bmatrix} = \vec 0 - all basic variables should be removed from the objective function. To do this, we'll take a linear combination \vec y of the equality constraints, to get {\vec y}^T \begin{bmatrix} 1 & 0 & 1 & -1 \\ 0 & 1 & 1 & 2 \end{bmatrix} \vec x = {\vec y}^T \begin{bmatrix} 1 \\ 2 \end{bmatrix}, or 0 = -{\vec y}^T \begin{bmatrix} 1 & 0 & 1 & -1 \\ 0 & 1 & 1 & 2 \end{bmatrix} \vec x + {\vec y}^T \begin{bmatrix} 1 \\ 2 \end{bmatrix}.
We'll add the objective function to both sides of that, to get \begin{bmatrix} 0 & 0 & 2 & 4 \end{bmatrix} \vec x + {\vec y}^T \begin{bmatrix} 1 & 0 & 1 & -1 \\ 0 & 1 & 1 & 2 \end{bmatrix} \vec x = \begin{bmatrix} 0 & 0 & 2 & 4 \end{bmatrix} \vec x + {\vec y}^T \begin{bmatrix} 1 \\ 2 \end{bmatrix}. Rearranging, we get \begin{bmatrix} 0 & 0 & 2 & 4 \end{bmatrix} \vec x = \left(\begin{bmatrix} 0 & 0 & 2 & 4 \end{bmatrix} - {\vec y}^T \begin{bmatrix} 1 & 0 & 1 & -1 \\ 0 & 1 & 1 & 2 \end{bmatrix}\right) \vec x + {\vec y}^T \begin{bmatrix} 1 \\ 2 \end{bmatrix}.
Clearly, the left side is the objective function, and it's equal to the right side. However, we can manipulate \vec y until the value that's being multiplied with \vec x has zero entries for our basic variables! Basically, we want to choose \vec y such that {\vec y}^T \begin{bmatrix} 1 & 0 & 1 & -1 \\ 0 & 1 & 1 & 2 \end{bmatrix}_{\set{2, 3}} = \begin{bmatrix} 0 & 0 & 2 & 4 \end{bmatrix}_{\set{2, 3}}. In other words, y_1 0 + y_2 1 = 0, y_1 1 + y_2 1 = 2. Solving this system, we get \vec y = \begin{bmatrix} 2 \\ 0 \end{bmatrix}.
Now \begin{bmatrix} 0 & 0 & 2 & 4 \end{bmatrix} \vec x = \left(\begin{bmatrix} 0 & 0 & 2 & 4 \end{bmatrix} - {\vec y}^T \begin{bmatrix} 1 & 0 & 1 & -1 \\ 0 & 1 & 1 & 2 \end{bmatrix}\right) \vec x + {\vec y}^T \begin{bmatrix} 1 \\ 2 \end{bmatrix} becomes \begin{bmatrix} 0 & 0 & 2 & 4 \end{bmatrix} \vec x = \left(\begin{bmatrix} 0 & 0 & 2 & 4 \end{bmatrix} - \begin{bmatrix} 2 & 0 & 2 & -2 \end{bmatrix}\right) \vec x + 2 = \begin{bmatrix} -2 & 0 & 0 & 6 \end{bmatrix} \vec x + 2. So for feasible solutions, the objective function is exactly equal to \begin{bmatrix} -2 & 0 & 0 & 6 \end{bmatrix} \vec x + 2.
In general, we take an arbitrary linear combination of the constraints {\vec y}^T A \vec x = {\vec y} \vec b, and then add the objective function to both sides to get {\vec c}^T \vec x + {\vec y}^T A \vec x = {\vec c}^T \vec x + {\vec y} \vec b, and rearrange to get {\vec c}^T \vec x = \left({\vec c}^T - {\vec y}^T A\right) \vec x + {\vec y} \vec b.
Then, we solve the linear system {\vec c}_B - {\vec y}^T A_B = {\vec 0}^T for \vec y. We can transpose both sides to write the formula as A_B^T \vec y = {\vec c}_B, then multiply both sides by (A_B^{-1})^T to get \vec y = (A_B^{-1})^T {\vec c}_B (this is equivalent because A_B is guaranteed to be non-singular). Then, we subsitute that value of \vec y back into {\vec c}^T \vec x = \left({\vec c}^T - {\vec y}^T A\right) \vec x + {\vec y} \vec b to get another way to write the objective function for feasible values of \vec x.
So in general, given the LP in SEF "maximize {\vec c}^T \vec x subject to A \vec x = \vec b and \vec x \ge \vec 0", the corresponding LP in canonical form for a basis B can be computed as follows:
Aside: sometimes in this course we'll use the notation A^{-T}, it represents inverse of transpose ((M^T)^{-1}) or transpose of inverse ((M^{-1})^T), which are equivalent to each other.
Suppose we have an LP "maximize {\vec c}^T \vec x subject to A \vec x = \vec b and \vec x \ge \vec 0" in canonical form for a basis B, and a basic feasible solution \vec x. The simplex algorithm requires us to find a better feasible solution from this, but how do we do that?
Essentially, to get a better feasible solution from an LP:
So, we started with an LP in canonical form for a feasible basis, and used it to obtain a better feasible basis (or a proof of unboundedness). We can then rewrite the LP in canonical form again and repeat this process until we either find an optimal solution or prove that it's unbounded. This is the simplex algorithm.
To summarize the procedure of finding an improved feasible basis given an existing feasible basis for an LP in canonical form:
This is essentially the entire simplex algorithm - it takes in an LP in SEF and a feasible basis, and outputs an optimal solution or proof that the LP is unbounded. To actually use this to solve arbitrary LPs though, we need to find a feasible basis first. This will be the next topic.
The only thing left to do is to ensure that the simplex algorithm terminates - how do we know it won't get stuck in a cycle of bases? To do this, we can use Bland's Rule, where whenever we need to choose a variable the leftmost variable in each case (the variable x_i with the smallest possible index i). This affects which variable we choose in step 4 and 7 above.
;wip: Why does Bland's rule work? proof outlined in exercise 10, page 69 of textbook
One of the interesting results we found is that the optimal solution given by the simplex algorithm is always a basic feasible solution, and that the simplex algorithm, with Bland's rule always gives us the optimal solution if there is one. LPs can have infinite feasible solutions, but the number of bases is finite, and so there are a finite number of basic feasible solutions.
This tells us that it's actually possible to find optimal solutions LPs by just enumerating every basis, finding the basic solution for each, and choosing the one that's both feasible and has the highest objective function value. However, this would be really inefficient - the number of possible bases goes up exponentially with respect to the number of constraints. That said, the simplex algorithm can also display exponential behaviour for LPs like the Klee-Minty cube, though it's a lot more efficient on real-world LPs.
One of the first steps before the simplex algorithm is, given an LP in SEF, to find a feasible basis, or detect that there is none. How do we find a basic feasible solution for an arbitrary LP in SEF? In other words, how do we find an \vec x \ge 0 such that A \vec x = \vec b (the objective function doesn't matter since we're just trying to find a basic feasible solution)?
Now we'll construct an auxilary problem - an LP that has a known feasible basis, whose solutions can be converted into a feasible basic solution for the original LP. We can then solve the auxilary problem using the simplex algorithm (since we have the auxilary LP and a feasible basis for that auxilary LP), and then convert the optimal solution to that auxilary LP into a basic feasible solution for our original LP.
The auxilary LP for an LP in SEF "maximize {\vec c}^T \vec x subject to A \vec x = \vec b and \vec x \ge \vec 0" (where A is an n-row-m-column matrix) can be constructed as follows:
Now we have the auxilary LP and a feasible basis for the auxilary LP B = \set{m + 1, \ldots, m + n}, so we can find an optimal solution using the simplex algorithm. Clearly, this LP is bounded, because it's a minimization problem where the objective function has all positive coefficients. Therefore, since we know that the LP isn't infeasible (since we have a feasible solution already), and the LP isn't unbounded, it must have an optimal solution!
The optimal solution to the auxilary LP gives us a feasible solution to the original LP, or tells us the original LP is infeasible. If any of the error variables x_{m + 1}, \ldots, x_{m + n} are non-zero, that means the original LP is infeasible, since if the original LP was feasible, it could've satisfied the constraints with all of the error variables set to 0, which would result in a smaller objective function value than we got.
Likewise, if all of the error variables are zero, then we satisfied the constraints without using any of the error variables, so \begin{bmatrix} x_1 \\ \vdots \\ x_m \end{bmatrix} is a feasible solution for the original LP. Additionally, since it's the output of the simplex algorithm, it's also a feasible basic solution, so we can construct a feasible basis from \set{i : x_i > 0, 1 \le i \le n}. This isn't guaranteed to give us n columns, but we can add in other column indices from the original LP until we do get a basis, checking them for linear independence each time.
Therefore, we can solve an arbitrary LP as follows:
This is called the two-phase algorithm.
Now we can actually prove the fundamental theorem of LPs using the simplex algorithm! Since we proved without using that theorem that the simplex algorithm works, and the simplex algorithm outputs either an optimal solution, proof of unboundedness, or proof of infeasibility, we've proved that an arbitrary LP must fall into one of those three categories.
In the geometric interpretation of LPs, the feasible region is the space of all possible solutions, \mb{R}^m where m is the number of variables in the LP.
A polyhedron P is a set of vectors that can be defined by P = \set{\vec x : A \vec x \le \vec b}, where A is a matrix and \vec b is a vector. Think of it as a solid n-dimensional shape - a bounded area of the n-dimensional space.
The feasible region of an LP is always a polyhedron. Given an LP "maximize {\vec c}^T \vec x subject to A \vec x = \vec b and \vec x \ge \vec 0", we can write its feasible region as \set{\vec x: \begin{bmatrix} A \\ -A \\ -I \end{bmatrix} \le \begin{bmatrix} \vec b \\ -\vec b \\ \vec 0 \end{bmatrix}}. Basically, this is just saying A \vec x \le \vec b, A \vec x \ge \vec b, and \vec x \ge 0.
A hyperplane is a set of vectors that can be defined by H = \set{\vec x : \vec a \cdot \vec x = \beta}, where \vec a is a non-zero vector and \beta \in \mb{R}. Think of it as an n-dimensional plane - an n - 1-dimensional space within the n-dimensional space. This is by definition the solution to a linear equation.
A halfspace is a set of vectors that can be defined by F = \set{\vec x : \vec a \cdot \vec x \le \beta}, much like a hyperplane. Think of it as the space on one side of a hyperplane, or half of the n-dimensional space. This is by definition the solution to a linear inequality.
A polyhedron can be thought of as an intersection of a finite number of halfspaces - each linear constraint can be thought of as a halfspace where points in the halfspace are those that satisfy that particular constraint, and points outside don't.
Recall that \vec x \cdot \vec y = \magn{x} \magn{y} \cos(\theta), where \theta is the minimum angle between \vec x and \vec y.
The hyperplane H = \set{\vec x : \vec a \cdot \vec x = 0} is the set of all vectors that are orthogonal to \vec a - the ones for which \theta = 90 \deg so \cos \theta = 0. In the same way, the halfspace F = \set{\vec x : \vec a \cdot \vec x \le 0} is the set of all vectors that are on the side of H that doesn't contain \vec a - the ones for which \theta \ge 90 \deg so \cos \theta \le 0.
The translate of a set of vectors S for some vector \vec p is the set \set{\vec x + \vec p : \vec x \in S}.
The hyperplane \set{\vec x : \vec a \cdot \vec x = \beta} is a translate of the hyperplane \set{\vec x : \vec a \cdot \vec x = 0}, by any vector \vec p such that \vec a \cdot \vec p = \beta.
The dimension of a hyperplane \set{\vec x : \vec a \cdot \vec x = \beta} is simply the dimension of \set{\vec x : \vec a \cdot \vec x = 0}, which is a vector space - its dimension is n - \mathrm{rank}(\vec a) = n - 1.
A line through two vectors \vec x, \vec y is defined as \set{\vec x + k(\vec y - \vec x), k \in \mb{R}}. A line segment is a subset of that, defined as \set{\vec x + k(\vec y - \vec x), 0 \le k \le 1}.
A set S \subseteq \mb{R}^n is convex if and only if, for all \vec x, \vec y \in S, the line segment between \vec x and \vec y is a subset of S (i.e., \forall x \in S, \forall y \in S, \set{\vec x + k(\vec y - \vec x), k \in \mb{R}} \subseteq S).
The union of two convex sets S_1 \cup S_2 is not necessarily convex - consider S_1 = \set{0}, S_2 = \set{1}, for example. However, the intersection of two convex sets S_1 \cap S_2 is convex. This is easily proved: any two points \vec x, \vec y in S_1 \cap S_2 is in both S_1 and S_2, by definition. Therefore, the line segment between them is a subset of both S_1 and S_2, and so the line segment is a subset of S_1 \cap S_2, as required.
All halfspaces are convex. This is also easy to prove:
Let F = \set{\vec x : \vec a \cdot \vec x \le \beta} be a halfspace.
Let \vec x_1, \vec x_2 \in F, and \vec v \in \set{\vec x_1 + k(\vec x_2 - \vec x_1), 0 \le k \le 1} be an arbitrary point on the line segment between \vec x_1 and \vec x_2.
So \vec a \cdot \vec v = \vec a \cdot \vec x_1 + \vec a \cdot k(\vec x_2 - \vec x_1) = (1 - k) \vec a \cdot \vec x_1 + k \vec a \cdot \vec x_2 for some 0 \le k \le 1.
Since \vec x_1 and \vec x_2 are in the halfspace, \vec a \cdot \vec x_2 \le \beta and \vec a \cdot \vec x_2 \le \beta \le \beta.
So \vec a \cdot \vec v = (1 - k) \vec a \cdot \vec x_1 + k \vec a \cdot \vec x_2 \le (1 - k) \beta + k \beta = \beta.
Therefore, \vec v \in F, and since \vec v is an arbitrary element of the line segment, the line segment must be a subset of F, which means that F is convex.
As mentioned earlier, a polyhedron is an intersection of a finite number of halfspaces. Since halfspaces are convex, and intersections of convex sets are convex, all polyhedra are convex. Since the feasible region of an LP is a polyhedron, it must also be convex - this is actually why solving LPs is relatively easy!
A point \vec v is properly contained in a line segment L = \set{\vec x + k(\vec y - \vec x), 0 \le k \le 1} if and only if \vec v \in L - \set{\vec x_1, \vec x_2} - when it's part of the line segment but isn't an endpoint.
Suppose we have a convex set S. A point \vec v \in S is an extreme point if and only if there is no line segment L \subseteq S that properly contains \vec v - if the point has to be one of the endpoints of any line segment in S that contains it.
Intuitively, an extreme point is one that's on a "curve" on the surface of the set. For a convex set that looks like a 2D square, the extreme points at the four corners. For a circle, every point is an extreme point - a convex set doesn't necessarily have a finite number of extreme points! A halfspace actually has no extreme points.
Suppose we have an LP with the constraints x_1 \le 3, x_2 \le 2 and x_1, x_2 \ge 0. The solution space is actually 2D for this LP, so each constraint forms a 2D halfspace (a hemiplane). If we plot the feasible region of this LP, we actually get a 3 by 2 square where the bottom left corner is the origin. Here, the extreme points are \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} 3, 0 \end{bmatrix}, \begin{bmatrix} 3 \\ 2 \end{bmatrix}, \begin{bmatrix} 2 \\ 0 \end{bmatrix}.
What's interesting is that all of these points, and only these points, satisfy two different constraints with equality at the same time. For example, for \begin{bmatrix} 3 \\ 2 \end{bmatrix}, the constraints x_1 \le 2 and x_2 \le 2 are satisfied with equality, so x_1 = 2 and x_2 = 2. In fact, this generalizes into multiple variables.
For an LP, given a constraint \vec p \cdot \vec x \le k and a solution \vec x, the constraint is tight if and only if it is satisfied with equality by \vec x - when \vec p \cdot x = k.
For an LP with constraints A \vec x \le \vec b, the tight constraints for a solution \vec x are denoted \overline A \vec x \le \overline{\vec b}, where \overline A and \overline{\vec b} are the rows of A and the entries of \vec b corresponding to constraints that are tight for \vec x.
Now we have the tools to define what an extreme point of a polyhedron is in a better way: given a point \vec x in a polyhedron \set{\vec x \in \mb{R}^m : A \vec x \le \vec b}, \vec x is an extreme point if and only if \mathrm{rank}(\overline A) = m. In other words, there's a theorem that says that an extreme point for an m-dimensional polyhedron is a point in the polyhedron that satisfies m linearly independent constraints with equality.
Now we'll prove the above:
Assume \mathrm{rank}(\overline A) = m. Suppose we have an arbitrary point \vec x in an arbitrary polyhedron P = \set{\vec x \in \mb{R}^m : A \vec x \le \vec b}.
Suppose \vec x is not an extreme point. Then \vec x is properly contained by a line segment L = \set{\vec x + k(\vec y - \vec x), 0 \le k \le 1} such that L \subseteq P, with some endpoints \vec v_1, \vec v_2. So there exists a 0 < k < 1 such that \vec x = \vec v_1 + k(\vec v_2 - \vec v_1) = k\vec v_2 + (1 - k)\vec v_1.
Clearly, the tight constraints \overline{\vec b} = \overline A \vec x = \overline A (k\vec v_2 + (1 - k)\vec v_1) = k \overline A \vec v_2 + (1 - k) \overline A \vec v_1.
Since L \subseteq P, \vec v_1 \in P and \vec v_2 \in P, \overline A \vec v_2 \le \overline{\vec b} and \overline A \vec v_1 \le \overline{\vec b}.
For any vectors \vec a, \vec b, \vec b, if \vec a = k \vec b + (1 - k) \vec c and \vec b \le \vec a and \vec c \le \vec a, then it must be that \vec a = \vec b = \vec c (if the weighted average of two numbers is greater or equal to both numbers, it must be equal to those numbers). This is easily proven by showing that \vec a = k\vec b + (1 - k)\vec c \le k\vec a + (1 - k)\vec a = \vec a, so k\vec b + (1 - k)\vec c = k\vec a + (1 - k)\vec a, so \vec b = \vec a and \vec c = \vec a.
Therefore, \overline{\vec b} = \overline A \vec v_1 = \overline A \vec v_2.
Since \mathrm{rank}(\overline A) = m, \overline A is invertible. Since \overline{\vec b} = \overline A \vec v_1 = \overline A \vec v_2, \vec v_1 = \vec v_2.
This is a contradiction, because then \vec x can't possibly be contained in L, if the line segment only has its endpoints. Therefore, if \mathrm{rank}(\overline A) = m, then \vec x must be an extreme point. Basically, we showed that if m independent constraints are satisfied by a point, then any line segment containing that point will only contain that point as both endpoints.
Assume \mathrm{rank}(\overline A) < m. So there must exist a vector \vec d such that \overline A \vec d = \vec 0.
Let \epsilon > 0. Let L be the line segment with endpoints \vec x - \epsilon \vec d and \vec x + \epsilon \vec d. Clearly, \vec x is properly contained by L, so we just need to show that L \subseteq P to prove that \vec x is not an extreme point.
Clearly, for the tight constraints \overline A, \overline A (\vec x - \epsilon \vec d) = \overline A \vec x - \epsilon \overline A \vec d. Since \overline A \vec d = \vec 0, \overline A \vec x - \epsilon \overline A \vec d = \overline A \vec x = \overline{\vec b}. Therefore, \vec x - \epsilon \vec d satisfies the tight constraints. The same argument shows that \vec x + \epsilon \vec d also satisfies the tight constraints.
Clearly, for any non-tight constraint \vec a \cdot \vec x \le \beta, \vec a \cdot (\vec x - \epsilon \vec d) = \vec a \cdot \vec x - \epsilon \vec a \cdot \vec d. Since \vec x \in P, we know that \vec a \cdot \vec x < \beta (the constraint is non-tight, so it isn't satisfied by equality, so we can use strict inequalities), but \vec a \cdot \vec d could be anything.
However, since \vec a \cdot \vec x < \beta, there must exist an \epsilon > 0 such that \vec a \cdot \vec x - \epsilon \vec a \cdot \vec d < \beta. Therefore, \vec x - \epsilon \vec d satisfies an arbitrary non-tight constraints, and by the same argument, so does \vec x + \epsilon \vec d.
Since all the tight constraints and arbitrary non-tight constraints are satisfied by \vec x - \epsilon \vec d and \vec x + \epsilon \vec d, so both are in P. Since P is convex and the endpoints are in P, L \subseteq P, so we have a line in P that properly contains \vec x.
Therefore, \vec x is not an extreme point if \mathrm{rank}(\overline A) < m.
Suppose we now have the feasible region of an LP P = \set{\vec x : \begin{bmatrix} 1 & 0 & -1 \\ 0 & 1 & 3 \end{bmatrix} \vec x = \begin{bmatrix} 2 \\ 4 \end{bmatrix}, \vec x \ge 0}. Recall that the feasible region of an LP is always a polyhedron. How do we check if \vec v = \begin{bmatrix} 2 \\ 4 \\ 0 \end{bmatrix} is an extreme point?
First, we can rewrite P in the form \set{\vec x : A \vec x = \vec b}: \set{\vec x : \begin{bmatrix} 1 & 0 & -1 \\ 0 & 1 & 3 \\ -1 & 0 & 1 \\ 0 & -1 & -3 \\ -1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & -1 \end{bmatrix} \vec x = \begin{bmatrix} 2 \\ 4 \\ -2 \\ -4 \\ 0 \\ 0 \\ 0 \end{bmatrix}}.
Now we take the constraints that are satified with equality by \vec v - namely, constraints 1, 2, 3, 4, and 7. So \overline A = \begin{bmatrix} 1 & 0 & -1 \\ 0 & 1 & 3 \\ -1 & 0 & 1 \\ 0 & -1 & -3 \\ 0 & 0 & -1 \end{bmatrix}. Clearly, \mathrm{rank}(\overline A) = 3, and the LP has 3 variables, so \vec v must be an extreme point.
An interesting consequence of the above theorem is that if we have the feasible region of an LP in SEF P = \set{\vec x : A \vec x = \vec b, \vec x \ge 0}, then \vec x \in P is an extreme point if and only if it's a basic feasible solution for that LP. In other words, the basic feasible solutions for an LP are the extreme points in its feasible region. This is easily proven from the previous theorem and the definition of a basic feasible solution.
Because the simplex algorithm moves between different bases until it reaches the optimal one, it's actually moving between extreme points in n-dimensional space. That means we could enumerate all the extreme points in the feasible region by enumerating all the basic feasible solutions.
Previously we looked at the shortest-path problem formulated as an IP - given a weighted graph G and s, t \in V(G), find the path from s to t that minimizes the total weight. Given a shortest-path problem and an s, t-path P, how do we prove if P is optimal? How do we find an optimal P, anyways?
The cardinality special case is the shortest path problem where all edge weights are 1 - we're trying to find the s, t-path P with the fewest edges. How do we show that P is an optimal solution, if it is?
Recall that any s, t-cut must contain one of the edges of P - if it didn't, then we could remove the edges in the cut and P would still connect s and t, which is a contradiction because removing the edges of an s, t-cut must disconnect s and t. As a result, if S \subseteq E(G) has an edge from every s, t-cut, then S contains an s, t-path.
If we can construct s, t-cuts \delta(U_1), \ldots, \delta(U_k) such that each one is disjoint from the other, then since an s, t-path must contain an edge from every s, t-cut, an s, t-cut must contain at least k edges. This allows us to lower bound the length of the shortest path, so if we find a path of length k, we know for sure that it is optimal. In other words, to prove that an s, t-path of length k is optimal, we find k pairwise disjoint s, t-cuts.
When we're looking at the general case (i.e., not the cardinality special case), we can do something similar. We construct a mapping between s, t-cuts and arbitrary non-negative values called cut widths: \overline y: \delta(U) \to \mb{R}_{\ge 0}. A feasible cut width assignment \overline y is such a mapping that, for every edge e, the sum of the widths of s, t-cuts containing e is less or equal to the weight of e. In other words, the cut widths summed up by edge shouldn't exceed the lengths of any edge. Formally, \forall e \in E(G), \sum_{\text{every s-t cut } U \text{ such that } e \in \delta(U)} \overline y_U \le c_e if and only if \overline y is feasible.
If we generalize the result of the cardinality special case, we find that if \overline y is a feasible cut width assignment, then every s, t-path must have length \sum_{\text{every s-t cut } U} \overline y_U or greater - the sum of the s, t-cut widths is a lower bound for the length of the shortest path! So if we find an , st path with that length, we know it's a shortest path.
Proof:
Let P be an s, t-path and \overline y be a feasible cut width assignment.
Clearly, the length of P is c(P) = \sum_{e \in P} c_e \ge \sum_{e \in P} \left(\sum_{\text{every s-t cut } \delta(U) \text{ such that } e \in \delta(U)} \overline y_U\right).
Since P is an s, t-path, P contains at least one edge from any \delta(U) (where U is an s, t-cut). So the condition for the inner sum will be satisfied, at one point or another, by every single s, t-cut.
So c(P) \ge \sum_{\text{every s-t cut } U} \overline y_U, giving us a lower bound for the length of P.
How do we find these cut width assignments, and do they always exist such that the total cut width is equal to the length of a shortest s, t-path? We will answer this using LPs.
Just as a feasible cut width assignment lower bounds the length of a shortest s, t-path, a feasible solution upper bounds the objective value of a minimization LP. Can we find a lower bound?
Suppose we have an LP "minimize {\vec c}^T \vec x subject to A \vec x \ge \vec b and \vec x \ge \vec 0". Clearly, any feasible solution \vec x must satisfy A \vec x \ge \vec b, so it must also satisfy {\vec k}^T A \vec x \ge {\vec k}^T \vec b for any \vec k \ge \vec 0 - we can construct any new constraint out of linear combinations of existing constraints, and \vec x will satisfy it.
Since {\vec c}^T \vec x \ge {\vec c}^T \vec x and (\vec k)^T A \vec x \ge (\vec k)^T \vec b, {\vec c}^T \vec x + (\vec k)^T A \vec x \ge {\vec c}^T \vec x + (\vec k)^T \vec b, so {\vec c}^T \vec x \ge (\vec k)^T \vec b + \left({\vec c}^T - (\vec k)^T A\right) \vec x.
Since we're trying to find a lower bound, we want to write the objective function value in the form z + \vec m \cdot \vec x, where \vec m \cdot \vec x \ge 0 - this means that z is a lower bound for the objective function value. To ensure (\vec k)^T \vec b + \left({\vec c}^T - (\vec k)^T A\right) \vec x is in this form, we want {\vec c}^T - (\vec k)^T A \ge \vec 0.
That means we're trying to choose \vec k such that {\vec c}^T \ge (\vec k)^T A. If we do find such a \vec k, then since \vec x \ge \vec 0, \left({\vec c}^T - (\vec k)^T A\right) \vec x \ge 0, so {\vec c}^T \vec x \ge {\vec c}^T \vec x \ge (\vec k)^T \vec b, and (\vec k)^T \vec b is a lower bound for the objective function.
We want the largest possible lower bound, to ensure that the actual minimum value isn't higher than our lower bound. To maximize this lower bound, we want to maximize (\vec k)^T \vec b subject to {\vec c}^T \ge (\vec k)^T A where \vec k \ge \vec 0, which is an LP such that the maximum objective function value is the minimum objective function value of the original LP. This new LP is the dual of the original LP. Since {\vec a}^T M = {M^T \vec a}^T for any vector \vec a and any matrix M, we can transpose both sides of the constraints to get "maximize {\vec b}^T \vec k subject to \vec c \ge A^T \vec k and \vec k \ge \vec 0".
In other words, "maximize {\vec b}^T \vec k subject to A^T \vec k \le \vec c and \vec k \ge \vec 0" is the dual LP for the primal LP "minimize {\vec c}^T \vec x subject to A \vec x \ge \vec b and \vec x \ge \vec 0", and vice versa. Basically, we transposed A, swapped \vec b and \vec c, and swapped the \ge and \le for all constraints except the non-negativity constraints.
Also, we can do the same thing to an LP in SEF to obtain its dual: "minimize {\vec b}^T \vec k subject to A^T \vec k \le \vec c" is the dual of "maximize {\vec c}^T \vec x subject to A \vec x = \vec b and \vec x \ge \vec 0".
The Weak Duality theorem says that if \vec x_1 is feasible for the primal maximization LP and \vec x_2 is feasible for the dual minimization LP, then \vec c \cdot \vec x_2 \le \vec b \cdot \vec x_1. This is because {\vec b}^T \vec x_2 = {\vec x_2}^T \vec b, and since \vec b \le A \vec x_1, {\vec b}^T \vec x_2 \le {\vec x_2}^T (A \vec x_1) = ({\vec x_2}^T A) \vec x_1 = (A^T \vec x_2)^T \vec x_1. Since A^T \vec x_2 \le \vec c, {\vec b}^T \vec x_2 \le {\vec c}^T \vec x_1, as required.
In other words, the objective function value for any feasible solution of the maximization LP is upper bounded by the objective function value for any feasible solution of the LP's dual, which is a minimizatoin problem. The main consequence of this is that if the objective function values are equal, both feasible solutions are optimal for their respective LPs.
As a result, the weak duality theorem also says that if \vec b \cdot \vec x_1 = \vec c \cdot \vec x_2 (the objective functions are exactly equal), then \vec x_1 is optimal for the primal LP and x_2 is optimal for the dual LP. This is because the dual's objective function bounds the primal's objective function, and vice versa, so if they're equal, they must be at their bounds.
The LP relaxation of an IP is just the IP with integer constraints removed.
Consider the shortest-path problem "minimize {\vec c}^T \vec x subject to A \vec x \ge \vec 1 and \vec x \ge \vec 0". The IP rows of A are s, t-cuts and columns are edges, and values of A are 1 or 0 depending on whether the corresponding edges are in the corresponding \delta(U) values for the s, t-cuts. Clearly, the solution \vec x for is a vector of integers 1 or 0, representing whether that given edge is a part of the shortest s, t-path, and we're trying to minimize the sum of weights of those edges.
Clearly, the objective function value is the length of the shortest s, t-path. Then the LP relaxation of the shortest-path IP must have an optimal value that's at most the optimal value for the IP, since removing constraints can only expand the feasible region.
Also note that the solution to the dual of the LP relaxation of the shortest-path IP is always a valid cut width assignment! This is because the dual, "maximize \vec 1 \cdot \vec y subject to A^T \vec y \le \vec c and \vec y \ge \vec 0", is maximizing the number of cuts while each constraint in A^T \vec y = \vec b ensures that cut weights summed for a given edge can't exceed the edge's weight. Each value of \vec y is a cut width, and each column of A^T is an s, t-cut. In other words, A[U, e] is 1 whenever e \in \delta(U) and 0 otherwise, where U is an s, t-cut and e is an edge.
From the Weak Duality theorem, the dual of the LP relaxation must have an optimal value that's less or equal to the optimal value for the LP relaxation. That means the cut width assignment \vec y has total width at most the length of the shortest s, t-path.
So if we find a feasible cut width assignment such that the total width is equal to the length of an s, t-path, then the s, t path must be a shortest s, t-path and the feasible cut width assignment must be maximal.
Also, the characteristic vector of an s, t-path P = sa, ab, bt is the vector \vec x where each element corresponds to an edge, and each element is 1 if the corresponding edge is in P, or 0 otherwise.
Week 8.
We looked at duals for the shortest path problem, but how did we find these shortest paths? How did we find the duals? In other words, how do we find an s, t-path P and a feasible dual (cut width assignment) such that the length of the path is equal to the total width of the cut width assignment?
In this course we will denote directed graphs as D = \tup{V, E}, and directed edges as \vec{uv}. A directed s, t-path might be denoted with \vec{sx}, \vec{xy}, \vec{yt}.
The slack of an edge given a feasible cut width assignment in a shortest-path problem is the difference between the edge's weight and the total width of the cuts that contain the edge. An equality edge is an edge with zero slack.
Also, an active s, t-cut in a feasible cut width assignment is an s, t-cut with a non-zero width.
Formally, the slack of an edge e for a feasible cut width assignment y is defined as \mathrm{slack}_y(e) = c_e - \sum_{\text{every st-cut } \delta(U) \text{ such that } e \in \delta(U)} y_U, where c_e is the weight of e and y_U is the width of the s, t-cut \delta(U).
There's actually a simple algorithm that will always find a shortest-path and a feasible cut width assignment of the same total width as the path's length (assuming non-negative edge weights):
Clearly, this algorithm must terminate, because every step adds a new vertex to U, so eventually it must reach t and terminate, assuming the graph is connected.
How do we know the resulting cut width assignment is of maximal width, or that the resulting s, t-path is of minimal length?
First of all, P is a shortest s, t-path if all edges in P are equality edges, and every active cut contains exactly one edge of P. Proof:
Assume all edges in P are equality edges, and every active cut contains exactly one edge of P.
Clearly, the length of the path is \sum_{e \in P} c_e. Since all edges are equality edges (so their slack is 0), c_e = \sum_{\text{every st-cut } \delta(U) \text{ such that } e \in \delta(U)} y_U for any edge e \in P.
So the length of the path is \sum_{e \in P} \sum_{\text{every st-cut } \delta(U) \text{ such that } e \in \delta(U)} y_U.
Clearly, if we expand the summation, e \in \delta(U) will occur \abs{P \cap \delta(U)} times for each s, t-cut \delta(U). So \sum_{e \in P} \sum_{\text{every st-cut } \delta(U) \text{ such that } e \in \delta(U)} y_U = \sum_{\text{every st-cut } \delta(U)} \abs{P \cap \delta(U)} y_U.
However, since every active s, t-cut contains exactly one edge of P, \abs{P \cap \delta(U)} = 1 for them, so \abs{P \cap \delta(U)} y_U = y_U. For any non-active s, t-cut U, y_U = 0, so \abs{P \cap \delta(U)} y_U = y_U. Therefore, \sum_{\text{every st-cut } \delta(U)} \abs{P \cap \delta(U)} y_U = \sum_{\text{every st-cut } \delta(U)} y_U.
Since \sum_{e \in P} c_e = \sum_{\text{every st-cut } \delta(U)} y_U, by Weak Duality, P is a minimum length path and y is a maximum width cut width assignment.
Now let's use this to prove that the algorithm gives the correct result. First, we note that at the end of every iteration in step 2, the following properties hold:
These invariants hold when the algorithm terminates. Therefore:
Since all edges in P are equality edges, and every active cut contains exactly one edge of P, we can use the proposition proved above to say that P is a shortest s, t-path. Therefore, the algorithm is always correct.
Note that if all edge weights are integers, then in the algorithm output, all cut widths are also all integers.
Week 9.
Previously, we looked at duality for the shortest-path problem. now we'll look at duality in general.
We can write an arbitrary LP as "maximize/minimize {\vec c}^T \vec x subject to A \vec x ? b and \vec x ? \vec 0". Here A \vec x ? b represents n individual constraints \mathrm{row}_i(A) \cdot \vec x \le b_i, \vec a_i \cdot \vec x \ge b_i, or \vec a_i \cdot \vec x = b_i. Likewise, \vec x ? \vec 0 represents n variable definitions x_i \ge 0, x_i \le 0, or x_i \text{ free}.
The dual has a variable for each constraint in the primal LP, and a constraint for each variable in the primal LP.
Recall that the dual of an LP "maximize {\vec c}^T \vec x subject to A \vec x \ge b and \vec x \ge \vec 0" is "minimize {\vec b}^T \vec y subject to A^T \vec y \le c and \vec y \ge \vec 0". This only works for LPs in a certain form, but it's easy to change LPs into this form:
Alternatively, a more mechanical way:
A good way to memorize maximization LP dual conversion is "convert constraints to variables (flipping inequalities in the process), convert variables to constraints, change maximize to minimize". A good way to memorize minimization LP dual conversion is "convert constraints to variables, convert variables to constraints (flipping inequalities in the process), change minimize to maximize".
We can generalize the weak duality theorem to these arbitrary LPs: if \vec x_1 is feasible for the primal maximization LP and \vec x_2 is feasible for the dual minimization LP, then \vec c \cdot \vec x_2 \le \vec b \cdot \vec x_1. Proof:
Clearly, all inequality constraints can be converted into equality constraints using slack variables (e.g., \vec a \cdot \vec x \ge b becomes \vec a \cdot \vec b + s = b with the slack variable s \le 0, and \vec a \cdot \vec x \le b becomes \vec a \cdot \vec b + s = b with the slack variable s \ge 0).
So all constraints in the primal LP can be written as A\vec x + \vec s = \vec b, with the slack variables represented with \vec s.
If we take the dual of this, we get the constraints A^t \vec y + \vec w = \vec c in the dual LP, where \vec w represents the slack variables.
Let \overline x be a feasible solution for the primal LP, and let \overline y be a feasible solution for the dual LP.
Let \overline s = \vec b - A \overline x represent the slack variable values in the primal LP, and let \overline w = \vec c - A^t \vec y represent the slack variable values in the dual LP.
Clearly, {\overline y}^T \vec b = {\overline y}^T (A \overline x + \overline s) = ({\overline y}^T A) \overline x + {\overline y}^T \overline s, since A \overline x = \overline b.
Clearly, ({\overline y}^T A) \overline x + {\overline y}^T \overline s = (\vec c - \overline w)^T \overline x + {\overline y}^T \overline s, since {\overline y}^T A = (A^T \overline y)^T and A^T \overline y + \overline w = \vec c, so {A^T \overline y}^T = (\vec c - \overline w)^T.
Clearly, (\vec c - \overline w)^T \overline x + {\overline y}^T \overline s = {\vec c}^T \overline x - {\overline w}^T \overline x + {\overline y}^T \overline s.
We know that {\overline w}^T \overline x \le 0 and {\overline y}^T \overline s \ge 0 (we won't prove this in this course, but it follows from the way we defined slack variables initially).
Therefore, {\overline y}^T \vec b \ge {\vec c}^T \overline x, as required.
This tells us the following:
The Strong Duality Theorem is similar to the weak duality, but it deals with optimal solutions rather than feasible solutions: if \vec x_1 is an optimal solution for a primal LP, then there exists a \vec x_2 that's an optimal solution for its dual LP, such that \vec c \cdot \vec x_2 = \vec b \cdot \vec x_1. Proof:
Assume our LP is in SEF, so it's of the form "maximize {\vec c}^T \vec x subject to A \vec x = \vec b and \vec x \ge \vec 0". Since any LP can be converted into SEF, this assumption is without loss of generality.
Then the dual must be "minimize {\vec b}^T \vec y subject to A^T \vec y \ge \vec c". Assume that the primal LP has an optimal solution \overline x.
Then the 2-phase simplex algorithm must terminate with an optimal basis B, so we can write the primal LP in canonical form for B and multiply both sides of the constraints equation by A^{-1}: "maximize z = {\overline y}^T \vec b - {\overline c}^T \vec x subject to \vec x_B = A_B^{-1} A_{\overline B} \vec x_{\overline B} and \vec x \ge \vec 0". Here, {\overline c}^T = {\vec c}^T - {\overline y}^T A.
So \overline x_B = A_B^{-1} \vec b and \overline x_{\overline B} = \vec 0. Note that the LP in canonical form is equivalent to the original LP in SEF - and \overline x has the same objective value in both.
Therefore, {\vec c}^T \overline x = {\overline y}^T \vec b = {\overline y}^T \vec b + {\vec c}^T \overline x.
Since \vec c has 0 for all of the basic components, {\vec c_B}^T \overline x_B = 0. Likewise, since \overline x_{\overline B} = \vec 0, {\vec c_{\overline B}}^T = 0. So {\vec c}^T \overline x = {\overline y}^T \vec b = \vec b \cdot \overline y.
Since B is an optimal basis for the , we know from the Simplex algorithm that \overline c \le \vec 0 in the canonical form objective function, so {\vec c}^T - {\overline y}^T A \le 0.
So {\vec c}^T \le {\overline y}^T A, which means \vec c \le A \overline y, which means that \overline y is feasible for the dual LP.
Since \overline y is a feasible solution whose objective function value is equal to that of a feasible solution of the primal LP, we've found an optimal solution to the dual LP that has the same objective function value, as required.
The strong duality theorem tells us that if a primal LP and its dual are both feasible, then they can neither be infeasible nor unbounded, so they must have optimal solutions that have equal value.
This gives us the slightly stronger version of the strong duality theorem: if a primal LP and its dual are both feasible, then both have optimal solutions with the same objective function value. Otherwise, one must be unbounded, and the other must be infeasible.
Recall the proof for the weak duality theorem, where we wrote out our LP as "maximize \vec c \cdot \vec x subject to A \vec x + \vec s = \vec b and \vec s \ge \vec 0", with the dual "minimize \vec b \cdot \vec y subject to A^T \vec y = \vec c and \vec y \ge 0". In the process, we showed that \overline y \cdot \vec b = \overline y \cdot (A \overline x + \overline s) = {\overline y}^T A \overline x + {\overline y}^T \overline s = (A^T \overline y)^T \overline x + {\overline y}^T \overline s = {\vec c}^T \overline x + {\overline y}^T \overline s.
According to the strong duality theorem, if \overline x, \overline y are optimal for their respective LPs, then {\vec c}^T \overline x = {\overline y}^T \vec b. So if \overline x, \overline y are optimal for their respective LPs and \overline y \cdot \vec b = {\vec c}^T \overline x + {\overline y}^T \overline s (\overline s represents the slack variables), then {\overline y}^T \overline s = 0.
Since \overline x, \overline y are feasible solutions, \overline x \ge 0 and \overline y \ge 0. Likewise, \overline s \ge 0 by definition. So if {\overline y}^T \overline s = 0, then either \overline y_i = 0 or \overline s_i = 0 for all 1 \le i \le n. In other words, either \overline y_i = 0, or constraint i in the primal LP is tight (satisfied with equality).
This gives us the complementary slackness special case: if \overline x is a feasible solution for a primal LP of the form "maximize \vec c \cdot \vec x subject to A \vec x \le \vec b" and \overline y is a feasible solution of its dual "minimize \vec b \cdot \vec y subject to A^T \vec y = \vec c and \vec y \ge 0", \overline x and \overline y are optimal for their respective LPs if and only if for each 1 \le i \le n, either \overline y_i = 0, or constraint i of the primal LP is tight.
The general complementary slackness theorem is the general version of the complementary slackness special case. It says that given a feasible solution \overline x of a primal LP and a feasible solution \overline y of its dual LP, \overline x and \overline y are optimal exactly when they satisfy the complementary slackness conditions:
Given \vec a_1, \ldots, \vec a_k \in \mb{R}^n, the cone generated by \vec a_1, \ldots, \vec a_k is C = \set{\lambda_1 \vec a_1 + \ldots + \lambda_k \vec a_k : \vec \lambda \ge \vec 0}. Essentially, it's the set of all points that can be generated by linear combinations of a set of vectors using non-negative coefficients.
The cone of tight constraints for an LP and a feasible solution \overline x is the cone generated by those rows of A that correspond to tight constraints for A \overline x = \vec b.
Interestingly, a feasible solution \overline x is an optimal solution if and only if \vec c is in the cone of tight constraints for \overline x. This follows from the general complementary slackness theorem, and is useful for geometrically interpreting what an optimal solution is.
Intuitively, we can think of \vec c as the "direction to optimize in". We're basically trying to move \overline x toward \vec c as far as possible. When \vec c is in the cone of tight constraints, \overline x must have been moved until it reached a polyhedron vertex, unable to move further past that vertex into the cone, which would make it optimal.
LPs are nice because we have ways of efficiently solving huge LP instances (millions of variables/constraints). Additionally, it's easy to construct short certificates for optimality (Strong Duality), feasibility (Farka's Lemma), and unboundedness. Plus, the fundamental theorem of LPs tells us there can only be three possible outcomes for any given LP.
In contrast, IPs don't have any known efficient solver algorithms, and we don't know of any way to generate short certificates for the solver outcomes. Plus, IPs have a lot more possible outcomes. On the other hand, IPs can be used to formulate many more types of practical problems, and certain special cases can easily be solved in practice.
For example, here's an IP that is feasible, bounded, and doesn't have an optimal solution: "maximize x_1 - \sqrt 2 x_2 subject to x_1 \le \sqrt 2 x_2, x_1 \ge 1, x_2 \ge 1, and x_1, x_2 \text{ integer}". This is because \sqrt 2 is irrational, so cannot be represented as a ratio. However, the IP tries to approximate \sqrt 2 as \frac{x_1}{x_2}, and the approximation can get arbitrarily better with larger integers x_1 and x_2. Therefore, no finite values for x_1, x_2 can be optimal. Formally:
Suppose x_1, x_2 is an optimal solution to the LP.
Let x_1' = 2x_1 + 2x_2, x_2' = x_1 + 2x_2. Clearly, x_1', x_2' satisfies constraint 2 if x_1, x_2 did.
Clearly, x_1' \le \sqrt 2 x_2' \iff 2x_1 + 2x_2 \le \sqrt 2 x_1 + \sqrt 2 2x_2 \iff (2 - \sqrt 2)x_1 \le \sqrt 2 2x_2 - 2 x_2 \iff (2 - \sqrt 2)x_1 \le (2 \sqrt 2 - 2) x_2 \iff x_1 \le \frac{2 \sqrt 2 - 2}{2 - \sqrt 2} x_2 = \sqrt 2 x_2.
Since x_1 \le \sqrt 2 x_2, x_1' \le \sqrt 2 x_2', so x_1', x_2' is a feasible solution.
Clearly, x_1' - \sqrt 2 x_2' > x_1 - \sqrt 2 x_2 \iff (2x_1 + 2x_2) - \sqrt 2(x_1 + 2x_2) > x_1 - \sqrt 2 x_2 \iff x_1 < \frac{2 - \sqrt 2}{\sqrt 2 - 1} x_2 = \sqrt 2 x_2.
Clearly, x_1 < \sqrt 2 x_2, because if x_1 = \sqrt 2 x_2, then \frac{x_1}{x_2} = \sqrt 2, which is impossible since x_1, x_2 are integers and \sqrt 2 is irrational. So x_1' - \sqrt 2 x_2' > x_1 - \sqrt 2 x_2.
So x_1, x_2 can't be optimal solutions to the LP, a contradiction. So the LP has no optimal solution.
In theory, IPs can be reduced to LPs, the LPs can be solved, and the solutions can be turned back into solutions for the IPs. Though this is procedure is itself no more efficient than existing ways to solve IPs, it tells us a lot about the relationship between IPs and LPs.
The convex hull of a subset C of \mb{R}^n is the smallest convex set that contains C. Intuitively, it's the shape formed by tightly wrapping the region in cling film - all gaps and concavities are filled in to form a convex region.
Also, it turns out that every subset C of \mb{R}^n has exactly one possible convex hull. This is because if we suppose there were two or more possible convex hulls for C, their set intersection is also a convex hull for C (because both sets being intersected are convex and contain C), but it would be smaller than at least one of the convex hulls, which contradicts the earlier fact that the two convex hulls are convex sets containing C of minimum size.
Meyer's Theorem: given the polyhedron P = \set{\vec x : A \vec x \le b} where all values in A and b are rational, the convex hull of all integer points in P is a polyhedron. The values must be rational because of cases like the one we saw above, where the IP's polyhedron is P = \set{\vec x : x_1 \le \sqrt 2 x_2, x_1 \ge 1, x_2 \ge 1}. Here, the convex hull of all integer points is not a polyhedron - if you try to construct one, it'll take infinitey many halfspaces.
Meyer's theorem means that solving IPs can be reduced to solving LPs. Suppose we have an IP "maximize \vec c \cdot \vec x subject to A \vec x \le \vec b and \vec x \text{ integer}", where all values in A, \vec b are rational. Clearly, the convex hull for the feasible solutions of the IP is a polyhedron A' \vec x \le \vec b', by Meyer's theorem.
Now consider the LP "maximize \vec c \cdot \vec x subject to A' \vec x \le \vec b'". This LP is infeasible exactly when the IP is infeasible, and is unbounded exactly when the IP is unbounded. Furthermore, an optimal solution to the IP is also an optimal solution to the LP, and an optimal solution to the LP that is also an extreme point is an optimal solution to the IP.
So to solve the IP, we could find A' and \vec b', use the Simplex algorithm to find an extreme optimal solution to the LP formed by A' and \vec b', and then use those as solutions as solutions to the IP. However, this is not a practical way of solving IPs, because the size of A' and \vec b' can grow exponentially with respect to the size of A and \vec b. One common real-world technique for solving IPs is to approximate the convex hull for the integer points.
First, we'll start by finding an optimal solution \overline x for the LP relaxation of the IP first (the LP relaxation is the IP with the integrality constraints removed). This is easy to do using the simplex algorithm, but the solution we get is not always an integral solution.
Now we'll add another constraint \vec \alpha \cdot \vec x \le \beta to the LP relaxation, such that the constraint is satisfied for all feasible solutions to the IP, but not for \overline x. Such a constraint is known as a cutting plane.
If we add a cutting plane constraint to our LP relaxation, and then find an optimal solution x', we note that the new optimal solution is now closer to at least one of the feasible solutions for the IP. By repeatedly adding cutting plane constraints and finding optimal solutions, we eventually run into an integer solution.
This process is called the cutting plane algorithm. Given the IP "maximize \vec c \cdot \vec x subject to A \vec x \le \vec b and x \text{ integer}":
It turns out the Simplex algorithm already finds cutting planes - take the final canonical form of the LP, then ;wip: 05b slides.
Recall that NLPs are problems of the form "maximize f(\vec x) subject to g_i(\vec x) \le 0 for all 1 \le i \le k". NLP feasible regions don't need to be convex, or even connected. This is what makes solving them hard.
We can also write this as "maximize \lambda subject to \lambda - f(\vec x) \le 0 and g_i(\vec x) \le 0 for all 1 \le i \le k". At an optimal solution, it is always true that \lambda = f(\vec x). Because \lambda is a linear function (of itself), we can always write any NLP such that it has a linear objective function. For the rest of this course, we'll assume that all NLPs we look at have linear objective functions.
NLPs generalize LPs and IPs. As we saw earlier, x_i = k can be written as x_i \le k and k \le x_i, x_i \text{ integer} can be written as \sin(\pi x_i) = 0, and 0 \le x_i \le 1, x_i \text{ integer} can be written as x_i(x_i - 1) = 0.
For our purposes, a feasible solution \vec x to a minimization problem is a local optimum if and only if there exists a \delta > 0 such that for any \vec x' such that \magn{\vec x' - \vec x} \le \delta, f(\vec x) < f(\vec x'). For a maximization problem, it's the same thing, but with f(\vec x) > f(\vec x') instead. Geometrically, local optima are feasible solutions for which there are no better solutions within some distance of it.
It turns out that if the feasible region is convex, all local optima are also global optima. Proof:
Let S be a convex feasible region for a minimization problem and \vec x \in S be local optimum. Suppose \vec x is not a global optimum. Let \vec x' \ne \vec x be a global optimum.
Since \vec x is a local optimum, for some \delta > 0, for any \vec z \in S such that \magn{\vec z - \vec x} \le \delta, \vec c \cdot \vec x \le \vec c \cdot \vec z.
Let \vec y = \lambda\vec x' + (1 - \lambda)\vec x. Since \vec x \in S and \vec x' \in S and S is convex, \vec y \in S.
Clearly, for \lambda = \delta, \magn{\vec y - \vec x} \le \delta. So \vec c \cdot \vec x \le \vec c \cdot \vec y.
However, \vec c \cdot \vec y = \vec c \cdot \left(\lambda\vec x' + (1 - \lambda)\vec x\right) = \lambda \vec c \cdot \vec x' + (1 - \lambda) \vec c \cdot \vec x.
Since \vec x' is a global optimum and \vec x isn't, \vec c \cdot \vec x' < \vec c \cdot \vec x. So \lambda \vec c \cdot \vec x' + (1 - \lambda) \vec c \cdot \vec x < \lambda \vec c \cdot \vec x + (1 - \lambda) \vec c \cdot \vec x = \vec c \cdot \vec x.
So, \vec c \cdot \vec y < \vec c \cdot \vec x and \vec c \cdot \vec x \le \vec c \cdot \vec y, a contradiction! Therefore, \vec x must be a global optimum.
Our strategy so far for finding globally optimal solutions to optimization problems has been to find a feasible solution, and then repeatedly refine it until we find one we know is locally optimal. For convex problems, locally optimal solutions are always globally optimal solutions, so this works well. However, since NLPs can have non-convex feasible regions, there can be locally optimal solutions that aren't globally optimal.
A convex function is a function f(\vec x) such that for all \vec a, \vec b, 0 \le k \le 1, f(k\vec a + (1 - k)\vec b) \le kf(\vec a) + (1 - k)f(\vec b). Basically, a convex function is one where the function's value at the interpolation of any two points is always less or equal to the interpolation of the function's value at those two points. For example, 1D convex functions like f(x) = x^2 are always flat or open upward. A common way of proving convexity is to show that a function's second derivative with respect to k is always positive.
For a convex function f(\vec x) and constant \beta, S = \set{\vec x \in \mb{R}^n : f(\vec x) \le \beta} is a convex set. Proof:
Let f(\vec x) be a convex function. Let \beta be a constant. Clearly, \vec x \in S if and only if f(\vec x) \le \beta.
Let \vec a, \vec b \in S, so f(\vec a) \le \beta and f(\vec b) \le \beta. Let \vec y = k\vec a + (1 - k)\vec b.
Since f(\vec x) is a convex function, f(\vec y) \le kf(\vec a) + (1 - k)f(\vec b). Clearly, kf(\vec a) + (1 - k)f(\vec b) \le k\beta + (1 - k)\beta = \beta.
Since \vec y \le \beta, S is a convex set.
The epigraph of a function f: \mb{R}^n \to \mb{R} is defined by \mathrm{epi}(f) = \set{\begin{bmatrix} y \\ \vec x \end{bmatrix} : y \ge f(\vec x), \vec x \in \mb{R}^n}. For 1D functions, it's the area above the curve; for 2D functions, it's the volume above the surface; and so on. The epigraph is a useful concept, because it's convex if and only if f(\vec x) is convex.
Suppose we have an NLP "maximize f(\vec x) subject to g_1(\vec x) \le 0 \land \ldots \land g_k(\vec x) \le 0". Let S_i = \set{\vec x \in \mb{R}^n : g_i(\vec x) \le 0} for all 1 \le i \le k. Clearly, the feasible region of this NLP is S_1 \cap \ldots \cap S_k - the intersection of all the values that satisfy each constraint.
If each g_i is convex, then so is each S_i, as we proved above. Since the intersection of convex sets is also convex, if each g_i is convex, so is the feasible region of the NLP.
The feasible region of an NLP is a region of \mb{R}^n, where n is the number of variables.
Suppose we have the NLP "minimize -x_1 - x_2 subject to -x_2 + x_1^2 \le 0 and -x_1 + x_2^2 \le 0 and -x_1 + \frac 1 2 \le 0". How do we prove that \overline x = \begin{bmatrix} 1 \\ 1 \end{bmatrix} is an optimal solution?
We can get a relaxation of the NLP by arbitrariy removing some constraints: "minimize -x_1 - x_2 subject to -x_2 + x_1^2 \le 0 and -x_1 + x_2^2 \le 0".
We can further relax it by replacing the -x_2 = x_1^2 \le 0 constraint with 2x_1 - x_2 \le 1, which is always satisfied when -x_2 = x_1^2 \le 0 is satisfied, and is satisfied by equality at \overline x.
We can do the same thing by replacing -x_1 + x_2^2 \le 0 with -x_1 + 2x_2 \le 1, which is also satisfied when -x_1 + x_2^2 \le 0 and satisfied by equality at \overline x.
The NLP then becomes "minimize -x_1 - x_2 subject to 2x_1 - x_2 \le 1 and -x_1 + 2x_2 \le 1", which is just an LP. Clearly, \overline x is an optimal solution to the LP relaxation, because the objective function coefficients are in the cone of tight constraints (we could also show this using strong duality, if we could find an optimal solution to the dual).
Since \overline x is an optimal solution to the LP relaxation, and is also a feasible solution to the NLP, it must also be an optimal solution to the NLP, since the LP's feasible region is a superset of the NLP's feasible region - each NLP feasible solution is an LP relaxation feasible solution, so the best of all the LP relaxation feasible solutions is also the best NLP feasible solution, if it is an NLP feasible solution.
We can apply this technique for showing that a feasible NLP solution is optimal to any NLP.
A subgradient of a convex function f: \mb{R}^n \to \mb{R} at a vector \overline x is a vector \vec s \in \mb{R}^n such that h(\vec x) \le f(\vec x) for all \vec x \in \mb{R}^n where h(\vec x) = f(\overline x) + \vec s \cdot (\vec x - \overline x).
Given a subgradient s, h(\vec x) is a linear function of \vec x, h(\overline x) = f(\overline x), and h(\vec x) is by definition always a lower bound for f(x). In other words, h(\vec x) is a linear lower bound for a convex function f(\vec x), and \vec s is the slope of h(\vec x).
Note that the subgradient is not unique - consider the subgradients possible at the extreme points of a polyhedron.
A halfspace F = \set{\vec x : \vec s \cdot \vec x \le \beta} is supporting a convex set C if and only if C \subseteq F (the halfspace contains the convex set) and \vec s \cdot \vec x = \beta (the halfspace is touching the convex set).
Given a convex function g(\vec x) and a vector \overline x such that g(\overline x) = 0 and a subgradient s of g at \overline x, F = \set{\vec x : h(\vec x) \le 0} is supporting C = \set{\vec x : g(\vec x) \le 0} where h(x) = g(\overline x) + \vec s \cdot (\vec x - \overline x).
Proof:
Let g(\vec x) be a convex function and \overline x be a value such that g(\overline x) = 0.
Let s be a subgradient of g at \overline x.
Since g(\vec x) is a convex function, C = \set{\vec x : g(\vec x) \le 0} is a convex set.
Since h(\vec x) is a linear function of \vec x, F = \set{\vec x : h(\vec x) \le 0} is a halfspace.
Since h(\overline x) = g(\overline x) = 0 = \beta, F is on the boundary of C.
Since h(\vec x) \le g(\vec x) \le 0, any \vec x \in F, so C \subseteq F.
So C is supporting F.
The above is useful because it allows us to construct LP relaxations of NLPs. Suppose we have the NLP "maximize f(\vec x) subject to g_1(\vec x) \le 0 \land \ldots \land g_k(\vec x) \le 0". If \overline x is a feasible solution of the NLP such that g_i(\overline x) = 0, g_i is convex, and \vec s is a subgradient of g_i at \overline x, then we can replace the constraint g_i(\vec x) \le 0 with h(\vec x) \le 0 to get a linear version of that constraint.
So given an NLP for which all g_i are convex, \overline x is a feasible solution such that g_i(\overline x) = 0 for all tight constraints i \in I, and \vec s_i is a subgradient for g_i(\vec x) at \overline x, \overline x is optimal -\vec c \in \mathrm{cone} \set{s_i : i \in I} (the negative objective function is in the cone of the tight constraints).
To get an LP relaxation of the NLP for a feasible solution \overline x, we can replace g_i(\vec x) with h_i(\vec x) for all tight constraints, and then discard the rest of the constraints.
Clearly, if the gradient \nabla f(\vec x) of a convex function f(\vec x) exists at \overline x, then it's a subgradient of f(\vec x) at \overline x. We compute the gradient of a function of a vector as \nabla f(\vec x) = \begin{bmatrix} \frac{\dee f(\vec x)}{\dee x_1} \\ \vdots \\ \frac{\dee f(\vec x)}{\dee x_n} \end{bmatrix}. Finding the gradient is a common way of finding a subgradient of a convex function at a given point.
A Slater point is a feasible solution \overline x of an NLP "minimize \vec c \cdot \vec x subject to g_i(\vec x) \le 0 for all 1 \le i \le k" such that for all 1 \le i \le k, g_i(\overline x) < 0 (a feasible solution that isn't on the boundary of the feasible region).
The Karush-Kuhn-Tucker theorem (KKT theorem) says that given an NLP "minimize \vec c \cdot \vec x subject to g_i(\vec x) \le 0 for all 1 \le i \le k" such that:
If this is the case, then \overline x is optimal if and only if -\vec c \in \mathrm{cone}\set{\nabla g_i(\overline x) : i \in I}.
The p-norm of a vector \vec x is defined as \magn{x}_p = \left(\sum \abs{x_i}^p\right)^{\frac 1 p}.
The infinity-norm of a vector \vec x is defined as \magn{x}_\infty = \max \set{\abs{x_1}, \ldots, \abs{x_n}}.
This work by Anthony Zhang is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
Copyright 2013-2017 Anthony Zhang.