|Making Hybrid Images||Neural Networks
|Bezier Curves and Picasso||Computing Homology||Probably Approximately
Correct – A Formal Theory
The standard inner product of two vectors has some nice geometric properties. Given two vectors , where by I mean the -th coordinate of , the standard inner product (which I will interchangeably call the dot product) is defined by the formula
This formula, simple as it is, produces a lot of interesting geometry. An important such property, one which is discussed in machine learning circles more than pure math, is that it is a very convenient decision rule.
In particular, say we’re in the Euclidean plane, and we have a line passing through the origin, with being a unit vector perpendicular to (“the normal” to the line).
If you take any vector , then the dot product is positive if is on the same side of as , and negative otherwise. The dot product is zero if and only if is exactly on the line , including when is the zero vector.
Here is an interactive demonstration of this property. Click the image below to go to the demo, and you can drag the vector arrowheads and see the decision rule change.
The code for this demo is available in a github repository.
It’s always curious, at first, that multiplying and summing produces such geometry. Why should this seemingly trivial arithmetic do anything useful at all?
The core fact that makes it work, however, is that the dot product tells you how one vector projects onto another. When I say “projecting” a vector onto another vector , I mean you take only the components of that point in the direction of . The demo shows what the result looks like using the red (or green) vector.
In two dimensions this is easy to see, as you can draw the triangle which has as the hypotenuse, with spanning one of the two legs of the triangle as follows:
If we call the (vector) leg of the triangle parallel to , while is the dotted line (as a vector, parallel to ), then as vectors . The projection of onto is just .
Another way to think of this is that the projection is , modified by removing any part of that is perpendicular to . Using some colorful language: you put your hands on either side of and , and then you squish onto along the line perpendicular to (i.e., along ).
And if is a unit vector, then the length of —that is, the length of the projection of onto —is exactly the inner product product .
Moreover, if the angle between and is larger than 90 degrees, the projected vector will point in the opposite direction of , so it’s really a “signed” length.
And this is precisely why the decision rule works. This 90-degree boundary is the line perpendicular to .
More technically said: Let be two vectors, and their dot product. Define by the length of , specifically . Define by by first letting , and then let . In words, you scale to a unit vector , use the result to compute the inner product, and then scale so that it’s length is . Then
Theorem: Geometrically, is the projection of onto the line spanned by .
This theorem is true for any -dimensional vector space, since if you have two vectors you can simply apply the reasoning for 2-dimensions to the 2-dimensional plane containing and . In that case, the decision boundary for a positive/negative output is the entire dimensional hyperplane perpendicular to (the projected vector).
In fact, the usual formula for the angle between two vectors, i.e. the formula , is a restatement of the projection theorem in terms of trigonometry. The part of the projection formula (how much you scale the output) is equal to . At the end of this post we have a proof of the cosine-angle formula above.
Part of why this decision rule property is so important is that this is a linear function, and linear functions can be optimized relatively easily. When I say that, I specifically mean that there are many known algorithms for optimizing linear functions, which don’t have obscene runtime or space requirements. This is a big reason why mathematicians and statisticians start the mathematical modeling process with linear functions. They’re inherently simpler.
In fact, there are many techniques in machine learning—a prominent one is the so-called Kernel Trick—that exist solely to take data that is not inherently linear in nature (cannot be fruitfully analyzed by linear methods) and transform it into a dataset that is. Using the Kernel Trick as an example to foreshadow some future posts on Support Vector Machines, the idea is to take data which cannot be separated by a line, and transform it (usually by adding new coordinates) so that it can. Then the decision rule, computed in the larger space, is just a dot product. Irene Papakonstantinou neatly demonstrates this with paper folding and scissors. The tradeoff is that the size of the ambient space increases, and it might increase so much that it makes computation intractable. Luckily, the Kernel Trick avoids this by remembering where the data came from, so that one can take advantage of the smaller space to compute what would be the inner product in the larger space.
Next time we’ll see how this decision rule shows up in an optimization problem: finding the “best” hyperplane that separates an input set of red and blue points into monochromatic regions (provided that is possible). Finding this separator is core subroutine of the Support Vector Machine technique, and therein lie interesting algorithms. After we see the core SVM algorithm, we’ll see how the Kernel Trick fits into the method to allow nonlinear decision boundaries.
Proof of the cosine angle formula
Theorem: The inner product is equal to , where is the angle between the two vectors.
Note that this angle is computed in the 2-dimensional subspace spanned by , viewed as a typical flat plane, and this is a 2-dimensional plane regardless of the dimension of .
Proof. If either or is zero, then both sides of the equation are zero and the theorem is trivial, so we may assume both are nonzero. Label a triangle with sides and the third side . Now the length of each side is and , respectively. Assume for the moment that is not 0 or 180 degrees, so that this triangle is not degenerate.
The law of cosines allows us to write
Moreover, The left hand side is the inner product of with itself, i.e. . We’ll expand using two facts. The first is trivial from the formula, that inner product is symmetric: . Second is that the inner product is linear in each input. In particular for the first input: and . The same holds for the second input by symmetry of the two inputs. Hence we can split up as follows.
Combining our two offset equations, we can subtract from each side and get
Which, after dividing by , proves the theorem if .
Now if or 180 degrees, the vectors are parallel, so we can write one as a scalar multiple of the other. Say for . In that case, . Now , since a norm is a length and is hence non-negative (but can be negative). Indeed, if are parallel but pointing in opposite directions, then , so , and . Otherwise and . This allows us to write , and this completes the final case of the theorem.
Problem: Determine if two polynomial expressions represent the same function. Specifically, if and are a polynomial with inputs, outputs and coefficients in a field , where is sufficiently large, then the problem is to determine if for every , in time polynomial in the number of bits required to write down and .
Solution: Let be the maximum degree of all terms in . Choose a finite set with . Repeat the following process 100 times:
- Choose inputs uniformly at random.
- Check if .
If every single time the two polynomials agree, accept the claim that they are equal. If they disagree on any input, reject. You will be wrong with probability at most .
Discussion: At first glance it’s unclear why this problem is hard.
If you have two representations of polynomials , say expressed in algebraic notation, why can’t you just do the algebra to convert them both into the same format, and see if they’re equal?
Unfortunately, that conversion can take exponential time. For example, suppose you have a polynomial . Though it only takes a few bits to write down, expressing it in a “canonical form,” often in the monomial form , would require exponentially many bits in the original representation. In general, it’s unknown how to algorithmically transform polynomials into a “canonical form” (so that they can be compared) in subexponential time.
Instead, the best we know how to do is treat the polynomials as black boxes and plug values into them.
Indeed, for single variable polynomials it’s well known that a nonzero degree polynomial has at most roots. A similar result is true for polynomials with many variables, and so we can apply that result to the polynomial to determine if . This theorem is so important (and easy to prove) that it deserves the name of lemma.
The Schwartz-Zippel lemma. Let be a nonzero polynomial of total degree over a field . Let be a finite subset of and let be chosen uniformly at random from . The probability that is at most .
Proof. By induction on the number of variables . For the case of , it’s the usual fact that a single-variable polynomial can have at most roots. Now for the inductive step, assume this is true for all polynomials with variables, and we will prove it for variables. Write as a polynomial in the variable , whose coefficients are other polynomials:
Here we’ve grouped by the powers of , so that are the coefficients of each . This is useful because we’ll be able to apply the inductive hypothesis to one of the ‘s, which have fewer variables.
Indeed, we claim there must be some which is nonzero for . Clearly, since is not the zero polynomial, some must be nonzero. If the only nonzero is , then we’re done because doesn’t depend on at all. Otherwise, take the largest nonzero . It’s true that the degree of is at most . This is true because the term has degree at most .
By the inductive hypothesis, if we choose and plug them into , we get zero with probability at most . The crucial part is that if this polynomial coefficient is nonzero, then the entire polynomial is nonzero. This is true even if an unlucky choice of causes the resulting evaluation .
To think about it a different way, imagine we’re evaluating the polynomial in phases. In the first phase, we pick the . We could also pick independently but not reveal what it is, for the sake of this story. Then we plug in the , and the result is a one-variable polynomial whose largest coefficient is . The inductive hypothesis tells us that this one-variable polynomial is the zero polynomial with probability at most . (It’s probably a smaller probability, since all the coefficients have to be zero, but we’re just considering the largest one for the sake of generality and simplicity)
Indeed, the resulting polynomial after we plug in has degree , so we can apply the inductive hypothesis to it as well, and the probability that it’s zero for a random choice of is at most .
Finally, the probability that both occur can be computed using basic probability algebra. Let be the event that, for these inputs, is zero, and the event that is zero for the and the additional .
Note the two quantities above that we don’t know are and , so we’ll bound them from above by 1. The rest of the quantities add up to exactly what we want, and so
which proves the theorem.
While this theorem is almost trivial to prove (it’s elementary induction, and the obvious kind), it can be used to solve polynomial identity testing, as well as finding perfect matchings in graphs and test numbers for primality.
But while the practical questions are largely solved–it’s hard to imagine a setting where you’d need faster primality testing than the existing randomized algorithms–the theory and philosophy of the result is much more interesting.
Indeed, checking two polynomials for equality has no known deterministic polynomial time algorithm. It’s one of a small class of problems, like integer factoring and the discrete logarithm, which are not known to be efficiently solvable in theory, but are also not known to be NP-hard, so there is still hope. The existence of this randomized algorithm increases hope (integer factorization sure doesn’t have one!). And more generally, the fact that there are so few natural problems in this class makes one wonder whether randomness is actually beneficial at all. From a polynomial-time-defined-as-efficient perspective, can every problem efficiently solvable with access to random bits also be solved without such access? In the computational complexity lingo, does P = BPP? Many experts think the answer is yes.
Problem: You have a catalog of items with discrete ratings (thumbs up/thumbs down, or 5-star ratings, etc.), and you want to display them in the “right” order.
Solution: In Python
''' score: [int], [int], [float] -> float Return the expected value of the rating for an item with known ratings specified by `ratings`, prior belief specified by `rating_prior`, and a utility function specified by `rating_utility`, assuming the ratings are a multinomial distribution and the prior belief is a Dirichlet distribution. ''' def score(self, ratings, rating_prior, rating_utility): ratings = [r + p for (r, p) in zip(ratings, rating_prior)] score = sum(r * u for (r, u) in zip(ratings, rating_utility)) return score / sum(ratings)
Discussion: This deceptively short solution can lead you on a long and winding path into the depths of statistics. I will do my best to give a short, clear version of the story.
As a working example I chose merely because I recently listened to a related podcast, say you’re selling mass-market romance novels—which, by all accounts, is a predictable genre. You have a list of books, each of which has been rated on a scale of 0-5 stars by some number of users. You want to display the top books first, so that time-constrained readers can experience the most titillating novels first, and newbies to the genre can get the best first time experience and be incentivized to buy more.
The setup required to arrive at the above code is the following, which I’ll phrase as a story.
Users’ feelings about a book, and subsequent votes, are independent draws from a known distribution (with unknown parameters). I will just call these distributions “discrete” distributions. So given a book and user, there is some unknown list of probabilities () for each possible rating a user could give for that book.
But how do users get these probabilities? In this story, the probabilities are the output of a randomized procedure that generates distributions. That modeling assumption is called a “Dirichlet prior,” with Dirichlet meaning it generates discrete distributions, and prior meaning it encodes domain-specific information (such as the fraction of 4-star ratings for a typical romance novel).
So the story is you have a book, and that book gets a Dirichlet distribution (unknown to us), and then when a user comes along they sample from the Dirichlet distribution to get a discrete distribution, which they then draw from to choose a rating. We observe the ratings, and we need to find the book’s underlying Dirichlet. We start by assigning it some default Dirichlet (the prior) and update that Dirichlet as we observe new ratings. Some other assumptions:
- Books are indistinguishable except in the parameters of their Dirichlet distribution.
- The parameters of a book’s Dirichlet distribution don’t change over time, and inherently reflect the book’s value.
So a Dirichlet distribution is a process that produces discrete distributions. For simplicity, in this post we will say a Dirichlet distribution is parameterized by a list of six integers , one for each possible star rating. These values represent our belief in the “typical” distribution of votes for a new book. We’ll discuss more about how to set the values later. Sampling a value (a book’s list of probabilities) from the Dirichlet distribution is not trivial, but we don’t need to do that for this program. Rather, we need to be able to interpret a fixed Dirichlet distribution, and update it given some observed votes.
The interpretation we use for a Dirichlet distribution is its expected value, which, recall, is the parameters of a discrete distribution. In particular if , then the expected value is a discrete distribution whose probabilities are
So you can think of each integer in the specification of a Dirichlet as “ghost ratings,” sometimes called pseudocounts, and we’re saying the probability is proportional to the count.
This is great, because if we knew the true Dirichlet distribution for a book, we could compute its ranking without a second thought. The ranking would simply be the expected star rating:
def simple_score(distribution): return sum(i * p for (i, p) in enumerate(distribution))
Putting books with the highest score on top would maximize the expected happiness of a user visiting the site, provided that happiness matches the user’s voting behavior, since the simple_score is just the expected vote.
Also note that all the rating system needs to make this work is that the rating options are linearly ordered. So a thumbs up/down (heaving bosom/flaccid member?) would work, too. We don’t need to know how happy it makes them to see a 5-star vs 4-star book. However, because as we’ll see next we have to approximate the distribution, and hence have uncertainty for scores of books with only a few ratings, it helps to incorporate numerical utility values (we’ll see this at the end).
Next, to update a given Dirichlet distribution with the results of some observed ratings, we have to dig a bit deeper into Bayes rule and the formulas for sampling from a Dirichlet distribution. Rather than do that, I’ll point you to this nice writeup by Jonathan Huang, where the core of the derivation is in Section 2.3 (page 4), and remark that the rule for updating for a new observation is to just add it to the existing counts.
Theorem: Given a Dirichlet distribution with parameters and a new observation of outcome , the updated Dirichlet distribution has parameters . That is, you just update the -th entry by adding to it.
This particular arithmetic to do the update is a mathematical consequence (derived in the link above) of the philosophical assumption that Bayes rule is how you should model your beliefs about uncertainty, coupled with the assumption that the Dirichlet process is how the users actually arrive at their votes.
The initial values for star ratings should be picked so that they represent the average rating distribution among all prior books, since this is used as the default voting distribution for a new, unknown book. If you have more information about whether a book is likely to be popular, you can use a different prior. For example, if JK Rowling wrote a Harry Potter Romance novel that was part of the canon, you could pretty much guarantee it would be popular, and set high compared to . Of course, if it were actually popular you could just wait for the good ratings to stream in, so tinkering with these values on a per-book basis might not help much. On the other hand, most books by unknown authors are bad, and should be close to zero. Selecting a prior dictates how influential ratings of new items are compared to ratings of items with many votes. The more pseudocounts you add to the prior, the less new votes count.
This gets us to the following code for star ratings.
def score(self, ratings, rating_prior): ratings = [r + p for (r, p) in zip(ratings, rating_prior)] score = sum(i * u for (i, u) in enumerate(ratings)) return score / sum(ratings)
The only thing missing from the solution at the beginning is the utilities. The utilities are useful for two reasons. First, because books with few ratings encode a lot of uncertainty, having an idea about how extreme a feeling is implied by a specific rating allows one to give better rankings of new books.
Second, for many services, such as taxi rides on Lyft, the default star rating tends to be a 5-star, and 4-star or lower mean something went wrong. For books, 3-4 stars is a default while 5-star means you were very happy.
The utilities parameter allows you to weight rating outcomes appropriately. So if you are in a Lyft-like scenario, you might specify utilities like [-10, -5, -3, -2, 1] to denote that a 4-star rating has the same negative impact as two 5-star ratings would positively contribute. On the other hand, for books the gap between 4-star and 5-star is much less than the gap between 3-star and 4-star. The utilities simply allow you to calibrate how the votes should be valued in comparison to each other, instead of using their literal star counts.
Hard to believe
Sanjeev Arora and his coauthors consider it “a basic tool [that should be] taught to all algorithms students together with divide-and-conquer, dynamic programming, and random sampling.” Christos Papadimitriou calls it “so hard to believe that it has been discovered five times and forgotten.” It has formed the basis of algorithms in machine learning, optimization, game theory, economics, biology, and more.
What mystical algorithm has such broad applications? Now that computer scientists have studied it in generality, it’s known as the Multiplicative Weights Update Algorithm (MWUA). Procedurally, the algorithm is simple. I can even describe the core idea in six lines of pseudocode. You start with a collection of objects, and each object has a weight.
Set all the object weights to be 1. For some large number of rounds: Pick an object at random proportionally to the weights Some event happens Increase the weight of the chosen object if it does well in the event Otherwise decrease the weight
The name “multiplicative weights” comes from how we implement the last step: if the weight of the chosen object at step is before the event, and represents how well the object did in the event, then we’ll update the weight according to the rule:
Think of this as increasing the weight by a small multiple of the object’s performance on a given round.
Here is a simple example of how it might be used. You have some money you want to invest, and you have a bunch of financial experts who are telling you what to invest in every day. So each day you pick an expert, and you follow their advice, and you either make a thousand dollars, or you lose a thousand dollars, or something in between. Then you repeat, and your goal is to figure out which expert is the most reliable.
This is how we use multiplicative weights: if we number the experts , we give each expert a weight which starts at 1. Then, each day we pick an expert at random (where experts with larger weights are more likely to be picked) and at the end of the day we have some gain or loss . Then we update the weight of the chosen expert by multiplying it by . Sometimes you have enough information to update the weights of experts you didn’t choose, too. The theoretical guarantees of the algorithm say we’ll find the best expert quickly (“quickly” will be concrete later).
In fact, let’s play a game where you, dear reader, get to decide the rewards for each expert and each day. I programmed the multiplicative weights algorithm to react according to your choices. Click the image below to go to the demo.
This core mechanism of updating weights can be interpreted in many ways, and that’s part of the reason it has sprouted up all over mathematics and computer science. Just a few examples of where this has led:
- In game theory, weights are the “belief” of a player about the strategy of an opponent. The most famous algorithm to use this is called Fictitious Play, and others include EXP3 for minimizing regret in the so-called “adversarial bandit learning” problem.
- In machine learning, weights are the difficulty of a specific training example, so that higher weights mean the learning algorithm has to “try harder” to accommodate that example. The first result I’m aware of for this is the Perceptron (and similar Winnow) algorithm for learning hyperplane separators. The most famous is the AdaBoost algorithm.
- Analogously, in optimization, the weights are the difficulty of a specific constraint, and this technique can be used to approximately solve linear and semidefinite programs. The approximation is because MWUA only provides a solution with some error.
- In mathematical biology, the weights represent the fitness of individual alleles, and filtering reproductive success based on this and updating weights for successful organisms produces a mechanism very much like evolution. With modifications, it also provides a mechanism through which to understand sex in the context of evolutionary biology.
- The TCP protocol, which basically defined the internet, uses additive and multiplicative weight updates (which are very similar in the analysis) to manage congestion.
- You can get easy -approximation algorithms for many NP-hard problems, such as set cover.
Additional, more technical examples can be found in this survey of Arora et al.
In the rest of this post, we’ll implement a generic Multiplicative Weights Update Algorithm, we’ll prove it’s main theoretical guarantees, and we’ll implement a linear program solver as an example of its applicability. As usual, all of the code used in the making of this post is available in a Github repository.
The generic MWUA algorithm
Let’s start by writing down pseudocode and an implementation for the MWUA algorithm in full generality.
In general we have some set of objects and some set of “event outcomes” which can be completely independent. If these sets are finite, we can write down a table whose rows are objects, whose columns are outcomes, and whose entry is the reward produced by object when the outcome is . We will also write this as for object and outcome . The only assumption we’ll make on the rewards is that the values are bounded by some small constant (by small I mean should not require exponentially many bits to write down as compared to the size of ). In symbols, . There are minor modifications you can make to the algorithm if you want negative rewards, but for simplicity we will leave that out. Note the table just exists for analysis, and the algorithm does not know its values. Moreover, while the values in are static, the choice of outcome for a given round may be nondeterministic.
The MWUA algorithm randomly chooses an object in every round, observing the outcome , and collecting the reward (or losing it as a penalty). The guarantee of the MWUA theorem is that the expected sum of rewards/penalties of MWUA is not much worse than if one had picked the best object (in hindsight) every single round.
Let’s describe the algorithm in notation first and build up pseudocode as we go. The input to the algorithm is the set of objects, a subroutine that observes an outcome, a black-box reward function, a learning rate parameter, and a number of rounds.
def MWUA(objects, observeOutcome, reward, learningRate, numRounds): ...
We define for object a nonnegative number we call a “weight.” The weights will change over time so we’ll also sub-script a weight with a round number , i.e. is the weight of object in round . Initially, all the weights are . Then MWUA continues in rounds. We start each round by drawing an example randomly with probability proportional to the weights. Then we observe the outcome for that round and the reward for that round.
# draw: [float] -> int # pick an index from the given list of floats proportionally # to the size of the entry (i.e. normalize to a probability # distribution and draw according to the probabilities). def draw(weights): choice = random.uniform(0, sum(weights)) choiceIndex = 0 for weight in weights: choice -= weight if choice <= 0: return choiceIndex choiceIndex += 1 # MWUA: the multiplicative weights update algorithm def MWUA(objects, observeOutcome, reward, learningRate numRounds): weights =  * len(objects) for t in numRounds: chosenObjectIndex = draw(weights) chosenObject = objects[chosenObjectIndex] outcome = observeOutcome(t, weights, chosenObject) thisRoundReward = reward(chosenObject, outcome) ...
Sampling objects in this way is the same as associating a distribution to each round, where if then the probability of drawing , which we denote , is . We don’t need to keep track of this distribution in the actual run of the algorithm, but it will help us with the mathematical analysis.
Next comes the weight update step. Let’s call our learning rate variable parameter . In round say we have object and outcome , then the reward is . We update the weight of the chosen object according to the formula:
In the more general event that you have rewards for all objects (if not, the reward-producing function can output zero), you would perform this weight update on all objects . This turns into the following Python snippet, where we hide the division by into the choice of learning rate:
# MWUA: the multiplicative weights update algorithm def MWUA(objects, observeOutcome, reward, learningRate, numRounds): weights =  * len(objects) for t in numRounds: chosenObjectIndex = draw(weights) chosenObject = objects[chosenObjectIndex] outcome = observeOutcome(t, weights, chosenObject) thisRoundReward = reward(chosenObject, outcome) for i in range(len(weights)): weights[i] *= (1 + learningRate * reward(objects[i], outcome))
One of the amazing things about this algorithm is that the outcomes and rewards could be chosen adaptively by an adversary who knows everything about the MWUA algorithm (except which random numbers the algorithm generates to make its choices). This means that the rewards in round can depend on the weights in that same round! We will exploit this when we solve linear programs later in this post.
But even in such an oppressive, exploitative environment, MWUA persists and achieves its guarantee. And now we can state that guarantee.
Theorem (from Arora et al): The cumulative reward of the MWUA algorithm is, up to constant multiplicative factors, at least the cumulative reward of the best object minus , where is the number of objects. (Exact formula at the end of the proof)
The core of the proof, which we’ll state as a lemma, uses one of the most elegant proof techniques in all of mathematics. It’s the idea of constructing a potential function, and tracking the change in that potential function over time. Such a proof usually has the mysterious script:
- Define potential function, in our case .
- State what seems like trivial facts about the potential function to write in terms of , and hence get general information about for some large .
- Theorem is proved.
- Wait, what?
Clearly, coming up with a useful potential function is a difficult and prized skill.
In this proof our potential function is the sum of the weights of the objects in a given round, . Now the lemma.
Lemma: Let be the bound on the size of the rewards, and a learning parameter. Recall that is the probability that MWUA draws object in round . Write the expected reward for MWUA for round as the following (using only the definition of expected value):
Then the claim of the lemma is:
Proof. Expand using the definition of the MWUA update:
Now distribute and split into two sums:
Using the fact that , we can replace with , which allows us to get
And then using the fact that (Taylor series), we can bound the last expression by , as desired.
Now using the lemma, we can get a hold on for a large , namely that
If then , simplifying the above. Moreover, the sum of the weights in round is certainly greater than any single weight, so that for every fixed object ,
Squeezing between these two inequalities and taking logarithms (to simplify the exponents) gives
Multiply through by , divide by , rearrange, and use the fact that when we have (Taylor series) to get
The bracketed term is the payoff of object , and MWUA’s payoff is at least a fraction of that minus the logarithmic term. The bound applies to any object , and hence to the best one. This proves the theorem.
Briefly discussing the bound itself, we see that the smaller the learning rate is, the closer you eventually get to the best object, but by contrast the more the subtracted quantity hurts you. If your target is an absolute error bound against the best performing object on average, you can do more algebra to determine how many rounds you need in terms of a fixed . The answer is roughly: let and pick . See this survey for more.
MWUA for linear programs
Now we’ll approximately solve a linear program using MWUA. Recall that a linear program is an optimization problem whose goal is to minimize (or maximize) a linear function of many variables. The objective to minimize is usually given as a dot product , where is a fixed vector and is a vector of non-negative variables the algorithm gets to choose. The choices for are also constrained by a set of linear inequalities, , where is a fixed vector and is a scalar for . This is usually summarized by putting all the in a matrix, in a vector, as
We can further simplify the constraints by assuming we know the optimal value in advance, by doing a binary search (more on this later). So, if we ignore the hard constraint , the “easy feasible region” of possible ‘s includes .
In order to fit linear programming into the MWUA framework we have to define two things.
- The objects: the set of linear inequalities .
- The rewards: the error of a constraint for a special input vector .
Number 2 is curious (why would we give a reward for error?) but it’s crucial and we’ll discuss it momentarily.
The special input depends on the weights in round (which is allowed, recall). Specifically, if the weights are , we ask for a vector in our “easy feasible region” which satisfies
For this post we call the implementation of procuring such a vector the “oracle,” since it can be seen as the black-box problem of, given a vector and a scalar and a convex region , finding a vector satisfying . This allows one to solve more complex optimization problems with the same technique, swapping in a new oracle as needed. Our choice of inputs, , are particular to the linear programming formulation.
Two remarks on this choice of inputs. First, the vector is a weighted average of the constraints in , and is a weighted average of the thresholds. So this this inequality is a “weighted average” inequality (specifically, a convex combination, since the weights are nonnegative). In particular, if no such exists, then the original linear program has no solution. Indeed, given a solution to the original linear program, each constraint, say , is unaffected by left-multiplication by .
Second, and more important to the conceptual understanding of this algorithm, the choice of rewards and the multiplicative updates ensure that easier constraints show up less prominently in the inequality by having smaller weights. That is, if we end up overly satisfying a constraint, we penalize that object for future rounds so we don’t waste our effort on it. The byproduct of MWUA—the weights—identify the hardest constraints to satisfy, and so in each round we can put a proportionate amount of effort into solving (one of) the hard constraints. This is why it makes sense to reward error; the error is a signal for where to improve, and by over-representing the hard constraints, we force MWUA’s attention on them.
At the end, our final output is an average of the produced in each round, i.e. . This vector satisfies all the constraints to a roughly equal degree. We will skip the proof that this vector does what we want, but see these notes for a simple proof. We’ll spend the rest of this post implementing the scheme outlined above.
Implementing the oracle
Fix the convex region for a known optimal value . Define as the problem of finding an such that .
For the case of this linear region , we can simply find the index which maximizes . If this value exceeds , we can return the vector with that value in the -th position and zeros elsewhere. Otherwise, the problem has no solution.
To prove the “no solution” part, say and you have a solution to . Then for whichever index makes bigger, say , you can increase without changing by replacing with and with zero. I.e., we’re moving the solution along the line until it reaches a vertex of the region bounded by and . This must happen when all entries but one are zero. This is the same reason why optimal solutions of (generic) linear programs occur at vertices of their feasible regions.
The code for this becomes quite simple. Note we use the
numpy library in the entire codebase to make linear algebra operations fast and simple to read.
def makeOracle(c, optimalValue): n = len(c) def oracle(weightedVector, weightedThreshold): def quantity(i): return weightedVector[i] * optimalValue / c[i] if c[i] > 0 else -1 biggest = max(range(n), key=quantity) if quantity(biggest) < weightedThreshold: raise InfeasibleException return numpy.array([optimalValue / c[i] if i == biggest else 0 for i in range(n)]) return oracle
Implementing the core solver
The core solver implements the discussion from previously, given the optimal value of the linear program as input. To avoid too many single-letter variable names, we use
linearObjective instead of .
def solveGivenOptimalValue(A, b, linearObjective, optimalValue, learningRate=0.1): m, n = A.shape # m equations, n variables oracle = makeOracle(linearObjective, optimalValue) def reward(i, specialVector): ... def observeOutcome(_, weights, __): ... numRounds = 1000 weights, cumulativeReward, outcomes = MWUA( range(m), observeOutcome, reward, learningRate, numRounds ) averageVector = sum(outcomes) / numRounds return averageVector
First we make the oracle, then the reward and outcome-producing functions, then we invoke the MWUA subroutine. Here are those two functions; they are closures because they need access to and . Note that neither nor the optimal value show up here.
def reward(i, specialVector): constraint = A[i] threshold = b[i] return threshold - numpy.dot(constraint, specialVector) def observeOutcome(_, weights, __): weights = numpy.array(weights) weightedVector = A.transpose().dot(weights) weightedThreshold = weights.dot(b) return oracle(weightedVector, weightedThreshold)
Implementing the binary search, and an example
Finally, the top-level routine. Note that the binary search for the optimal value is sophisticated (though it could be more sophisticated). It takes a max range for the search, and invokes the optimization subroutine, moving the upper bound down if the linear program is feasible and moving the lower bound up otherwise.
def solve(A, b, linearObjective, maxRange=1000): optRange = [0, maxRange] while optRange - optRange > 1e-8: proposedOpt = sum(optRange) / 2 print("Attempting to solve with proposedOpt=%G" % proposedOpt) # Because the binary search starts so high, it results in extreme # reward values that must be tempered by a slow learning rate. Exercise # to the reader: determine absolute bounds for the rewards, and set # this learning rate in a more principled fashion. learningRate = 1 / max(2 * proposedOpt * c for c in linearObjective) learningRate = min(learningRate, 0.1) try: result = solveGivenOptimalValue(A, b, linearObjective, proposedOpt, learningRate) optRange = proposedOpt except InfeasibleException: optRange = proposedOpt return result
Finally, a simple example:
A = numpy.array([[1, 2, 3], [0, 4, 2]]) b = numpy.array([5, 6]) c = numpy.array([1, 2, 1]) x = solve(A, b, c) print(x) print(c.dot(x)) print(A.dot(x) - b)
Attempting to solve with proposedOpt=500 Attempting to solve with proposedOpt=250 Attempting to solve with proposedOpt=125 Attempting to solve with proposedOpt=62.5 Attempting to solve with proposedOpt=31.25 Attempting to solve with proposedOpt=15.625 Attempting to solve with proposedOpt=7.8125 Attempting to solve with proposedOpt=3.90625 Attempting to solve with proposedOpt=1.95312 Attempting to solve with proposedOpt=2.92969 Attempting to solve with proposedOpt=3.41797 Attempting to solve with proposedOpt=3.17383 Attempting to solve with proposedOpt=3.05176 Attempting to solve with proposedOpt=2.99072 Attempting to solve with proposedOpt=3.02124 Attempting to solve with proposedOpt=3.00598 Attempting to solve with proposedOpt=2.99835 Attempting to solve with proposedOpt=3.00217 Attempting to solve with proposedOpt=3.00026 Attempting to solve with proposedOpt=2.99931 Attempting to solve with proposedOpt=2.99978 Attempting to solve with proposedOpt=3.00002 Attempting to solve with proposedOpt=2.9999 Attempting to solve with proposedOpt=2.99996 Attempting to solve with proposedOpt=2.99999 Attempting to solve with proposedOpt=3.00001 Attempting to solve with proposedOpt=3 Attempting to solve with proposedOpt=3 # note %G rounds the printed values Attempting to solve with proposedOpt=3 Attempting to solve with proposedOpt=3 Attempting to solve with proposedOpt=3 Attempting to solve with proposedOpt=3 Attempting to solve with proposedOpt=3 Attempting to solve with proposedOpt=3 Attempting to solve with proposedOpt=3 Attempting to solve with proposedOpt=3 Attempting to solve with proposedOpt=3 [ 0. 0.987 1.026] 3.00000000425 [ 5.20000072e-02 8.49831849e-09]
So there we have it. A fiendishly clever use of multiplicative weights for solving linear programs.
One of the nice aspects of MWUA is it’s completely transparent. If you want to know why a decision was made, you can simply look at the weights and look at the history of rewards of the objects. There’s also a clear interpretation of what is being optimized, as the potential function used in the proof is a measure of both quality and adaptability to change. The latter is why MWUA succeeds even in adversarial settings, and why it makes sense to think about MWUA in the context of evolutionary biology.
This even makes one imagine new problems that traditional algorithms cannot solve, but which MWUA handles with grace. For example, imagine trying to solve an “online” linear program in which over time a constraint can change. MWUA can adapt to maintain its approximate solution.
The linear programming technique is known in the literature as the Plotkin-Shmoys-Tardos framework for covering and packing problems. The same ideas extend to other convex optimization problems, including semidefinite programming.
If you’ve been reading this entire post screaming “This is just gradient descent!” Then you’re right and wrong. It bears a striking resemblance to gradient descent (see this document for details about how special cases of MWUA are gradient descent by another name), but the adaptivity for the rewards makes MWUA different.
Even though so many people have been advocating for MWUA over the past decade, it’s surprising that it doesn’t show up in the general math/CS discourse on the internet or even in many algorithms courses. The Arora survey I referenced is from 2005 and the linear programming technique I demoed is originally from 1991! I took algorithms classes wherever I could, starting undergraduate in 2007, and I didn’t even hear a whisper of this technique until midway through my PhD in theoretical CS (I did, however, study fictitious play in a game theory class). I don’t have an explanation for why this is the case, except maybe that it takes more than 20 years for techniques to make it to the classroom. At the very least, this is one good reason to go to graduate school. You learn the things (and where to look for the things) which haven’t made it to classrooms yet.
Until next time!