Making Hybrid Images

monalisa

The Mona Lisa

Leonardo da Vinci’s Mona Lisa is one of the most famous paintings of all time. And there has always been a discussion around her enigmatic smile. He used a trademark Renaissance technique called sfumato, which involves many thin layers of glaze mixed with subtle pigments. The striking result is that when you look directly at Mona Lisa’s smile, it seems to disappear. But when you look at the background your peripherals see a smiling face.

One could spend decades studying the works of these masters from various perspectives, but if we want to hone in on the disappearing nature of that smile, mathematics can provide valuable insights. Indeed, though he may not have known the relationship between his work and da Vinci’s, hundreds of years later Salvador Dali did the artist’s equivalent of mathematically isolating the problem with his painting, “Gala Contemplating the Mediterranean Sea.”

gala-dali

Gala Contemplating the Mediterranean Sea (Salvador Dali, 1976)

Here you see a woman in the foreground, but step back quite far from the picture and there is a (more or less) clear image of Abraham Lincoln. Here the question of gaze is the blaring focus of the work. Now of course Dali and da Vinci weren’t scribbling down equations and computing integrals; their artistic expression was much less well-defined. But we the artistically challenged have tools of our own: mathematics, science, and programming.

In 2006 Aude Oliva, Antonio Torralba, and Philippe. G. Schyns used those tools to merge the distance of Dali and the faded smiles of da Vinci into one cohesive idea. In their 2006 paper they presented the notion of a “hybrid image,” presented below.

monalisas

The Mona Lisas of Science

If you look closely, you’ll see three women, each of which looks the teensiest bit strange, like they might be trying to suppress a smile, but none of them are smiling. Blur your eyes or step back a few meters, and they clearly look happy. The effect is quite dramatic. At the risk of being overly dramatic, these three women are literally modern day versions of Mona Lisa, the “Mona Lisas of Science,” if you will.

Another, perhaps more famous version of their technique, since it was more widely publicized, is their “Marilyn Einstein,” which up close is Albert Einstein and from far away is Marilyn Monroe.

marilyn-einstein

Marilyn Einstein

This one gets to the heart of the question of what the eye sees at close range versus long range. And it turns out that you can address this question (and create brilliant works of art like the ones above) with some basic Fourier analysis.

Intuitive Fourier analysis (and references)

The basic idea of Fourier analysis is the idea that smooth functions are hard to understand, and realization of how great it would be if we could decompose them into simpler pieces. Decomposing complex things into simpler parts is one of the main tools in all of mathematics, and Fourier analysis is one of the clearest examples of its application.

In particular, the things we care about are functions f(x) with specific properties I won’t detail here like “smoothness” and “finiteness.” And the building blocks are the complex exponential functions

\displaystyle e^{2 \pi i kx}

where k can be any integer. If you have done some linear algebra (and ignore this if you haven’t), then I can summarize the idea succinctly by saying the complex exponentials form an orthonormal basis for the vector space of square-integrable functions.

Back in colloquial language, what the Fourier theorem says is that any function of the kind we care about can be broken down into (perhaps infinitely many) pieces of this form called Fourier coefficients (I’m abusing the word “coefficient” here). The way it’s breaking down is also pleasingly simple: it’s a linear combination. Informally that means you’re just adding up all the complex exponentials with specific weights for each one. Mathematically, the conversion from the function to its Fourier coefficients is called the Fourier transform, and the set of all Fourier coefficients together is called the Fourier spectrum. So if you want to learn about your function f, or more importantly modify it in some way, you can inspect and modify its spectrum instead. The reason this is useful is that Fourier coefficients have very natural interpretations in sound and images, as we’ll see for the latter.

We wrote f(x) and the complex exponential as a function of one real variable, but you can do the same thing for two variables (or a hundred!). And, if you’re willing to do some abusing and ignore the complexness of complex numbers, then you can visualize “complex exponentials in two variables” as images of stripes whose orientation and thickness correspond to two parameters (i.e., the k in the offset equation becomes two coefficients). The video below shows how such complex exponentials can be used to build up an image of striking detail. The left frame shows which complex exponential is currently being added, and the right frame shows the layers all put together. I think the result is quite beautiful.

This just goes to show how powerful da Vinci’s idea of fine layering is: it’s as powerful as possible because it can create any image! 

Now for digital images like the one above, everything is finite. So rather than have an infinitely precise function and a corresponding infinite set of Fourier coefficients, you get a finite list of sampled values (pixels) and a corresponding grid of Fourier coefficients. But the important and beautiful theorem is, and I want to emphasize how groundbreakingly important this is:

If you give me an image (or any function!) I can compute the decomposition very efficiently.

And the same theorem lets you go the other way: if you give me the decomposition, I can compute the original function’s samples quite easily. The algorithm to do this is called the Fast Fourier transform, and if any piece of mathematics or computer science has a legitimate claim to changing the world, it’s the Fast Fourier transform. It’s hard to pinpoint specific applications, because the transform is so ubiquitous across science and engineering, but we definitely would not have cell phones, satellites, internet, or electronics anywhere near as small as we do without the Fourier transform and the ability to compute it quickly.

Constructing hybrid images is one particularly nice example of manipulating the Fourier spectrum of two images, and then combining them back into a single image. That’s what we’ll do now.

As a side note, by the nature of brevity, the discussion above is a big disservice to the mathematics involved. I summarized and abused in ways that mathematicians would object to. If you want to see a much better treatment of the material, this blog has a long series of posts developing Fourier transforms and their discrete analogues from scratch. See our four primers, which lead into the main content posts where we implement the Fast Fourier transform in Python and use it to apply digital watermarks to an image. Note that in those posts, as in this one, all of the materials and code used are posted on this blog’s Github page.

High and low frequencies

For images, interpreting ranges of Fourier coefficients is easy to do. You can imagine the coefficients lying on a grid in the plane like so:

sherlock-spectrum

Each dot in this grid corresponds to how “intense” the Fourier coefficient is. That is, it’s the magnitude of the (complex) coefficient of the corresponding complex exponential. Now the points that are closer to the origin correspond informally to the broad, smooth changes in the image. These are called “low frequency” coefficients. And points that are further away correspond to sharp changes and edges, and are likewise called “high frequency” components. So the if you wanted to “hybridize” two images, you’d pick ones with complementary intensities in these regions. That’s why Einstein (with all his wiry hair and wrinkles) and Monroe (with smooth features) are such good candidates. That’s also why, when we layered the Fourier components one by one in the video from earlier, we see the fuzzy shapes emerge before the fine details.

Moreover, we can “extract” the high frequency Fourier components by simply removing the low frequency ones. It’s a bit more complicated than that, since you want the transition from “something” to “nothing” to be smooth in sone sense. A proper discussion of this would go into sampling and the Nyquist frequency, but that’s beyond the scope of this post. Rather, we’ll just define a family of “filtering functions” without motivation and observe that they work well.

Definition: The Gaussian filter function with variance \sigma and center (a, b) is the function

\displaystyle g(x,y) = e^{-\frac{(x - a)^2 + (y - b)^2}{2 \sigma^2}}

It looks like this

image credit Wikipedia

image credit Wikipedia

In particular, at zero the function is 1 and it gradually drops to zero as you get farther away. The parameter \sigma controls the rate at which it vanishes, and in the picture above the center is set to (0,0).

Now what we’ll do is take our image, compute its spectrum, and multiply coordinatewise with a certain Gaussian function. If we’re trying to get rid of high-frequency components (called a “low-pass filter” because it lets the low frequencies through), we can just multiply the Fourier coefficients directly by the filter values g(x,y), and if we’re doing a “high-pass filter” we multiply by 1 - g(x,y).

Before we get to the code, here’s an example of a low-pass filter. First, take this image of Marilyn Monroe

marilyn

Now compute its Fourier transform

dft

Apply the low-pass filter

filtered-dft

And reverse the Fourier transform to get an image

low-passed-marilyn

In fact, this is a common operation in programs like photoshop for blurring an image (it’s called a Gaussian blur for obvious reasons). Here’s the python code to do this. You can download it along with all of the other resources used in making this post on this blog’s Github page.

import numpy
from numpy.fft import fft2, ifft2, fftshift, ifftshift
from scipy import misc
from scipy import ndimage
import math

def makeGaussianFilter(numRows, numCols, sigma, highPass=True):
   centerI = int(numRows/2) + 1 if numRows % 2 == 1 else int(numRows/2)
   centerJ = int(numCols/2) + 1 if numCols % 2 == 1 else int(numCols/2)

   def gaussian(i,j):
      coefficient = math.exp(-1.0 * ((i - centerI)**2 + (j - centerJ)**2) / (2 * sigma**2))
      return 1 - coefficient if highPass else coefficient

   return numpy.array([[gaussian(i,j) for j in range(numCols)] for i in range(numRows)])

def filterDFT(imageMatrix, filterMatrix):
   shiftedDFT = fftshift(fft2(imageMatrix))
   filteredDFT = shiftedDFT * filterMatrix
   return ifft2(ifftshift(filteredDFT))

def lowPass(imageMatrix, sigma):
   n,m = imageMatrix.shape
   return filterDFT(imageMatrix, makeGaussianFilter(n, m, sigma, highPass=False))

def highPass(imageMatrix, sigma):
   n,m = imageMatrix.shape
   return filterDFT(imageMatrix, makeGaussianFilter(n, m, sigma, highPass=True))

if __name__ == "__main__":
   marilyn = ndimage.imread("marilyn.png", flatten=True)
   lowPassedMarilyn = lowPass(marilyn, 20)
   misc.imsave("low-passed-marilyn.png", numpy.real(lowPassedMarilyn))

The first function samples the values from a Gaussian function with the specified parameters, discretizing the function and storing the values in a matrix. Then the filterDFT function applies the filter by doing coordinatewise multiplication (note these are all numpy arrays). We can do the same thing with a high-pass filter, producing the edgy image below

high-passed-marilyn

And if we compute the average of these two images, we basically get back to the original.

sum-of-marilyns

So the only difference between this and a hybrid image is that you take the low-passed part of one image and the high-passed part of another. Then the art is in balancing the parameters so as to make the averaged image look right. Indeed, with the following picture of Einstein and the above shot of Monroe, we can get a pretty good recreation of the Oliva-Torralba-Schyns piece. I think with more tinkering it could be even better (I did barely any centering/aligning/resizing to the original images).

Albert Einstein, Marilyn Monroe, and their hybridization.

Albert Einstein, Marilyn Monroe, and their hybridization.

And here’s the code for it

def hybridImage(highFreqImg, lowFreqImg, sigmaHigh, sigmaLow):
   highPassed = highPass(highFreqImg, sigmaHigh)
   lowPassed = lowPass(lowFreqImg, sigmaLow)

   return highPassed + lowPassed

Interestingly enough, doing it in reverse doesn’t give quite as pleasing results, but it still technically works. So there’s something particularly important that the high-passed image does have a lot of high-frequency components, and vice versa for the low pass.

backwards

You can see some of the other hybrid images Oliva et al constructed over at their web gallery.

Next Steps

How can we take this idea further? There are a few avenues I can think of. The most obvious one would be to see how this extends to video. Could one come up with generic parameters so that when two videos are hybridized (frame by frame, using this technique) it is only easy to see one at close distance? Or else, could we apply a three-dimensional transform to a video and modify that in some principled way? I think one would not likely find anything astounding, but who knows?

Second would be to look at the many other transforms we have at our disposal. How does manipulating the spectra of these transforms affect the original image, and can you make images that are hybridized in senses other than this one?

And finally, can we bring this idea down in dimension to work with one-dimensional signals? In particular, can we hybridize music? It could usher in a new generation of mashup songs that sound different depending on whether you wear earmuffs :)

Until next time!

About these ads

AMS Network Science Mathematical Research Community

I don’t usually write promotional posts because I don’t enjoy reading them as much as I enjoy reading the technical posts. But I know that a lot of early graduate students and undergraduates read my blog, and this would be of interest to many of them.

I just got back from Utah yesterday where I attended a 5-day workshop run by the American Mathematical Society, called the Network Science Mathematical Research Community (MRC).

The point of the program is to bring graduate students and early career folks together from all over the country to start new collaborations. The AMS runs multiple MRC sessions every year, and this year the topics ranged from network science to quantum physics. We had a group of about 20 people, including statisticians, applied mathematicians, computer scientists, and a handful of pure combinatorialists. We self-organized into groups of four, and spent pretty much all day for the next four days eating great food, thinking about problems, proving theorems, enjoying the view, and discussing our ideas with the three extremely smart, successful, and amicable organizers. There were also career panels every evening that were, in my opinion, better than the average career panel.

The network science group (you can see me peeking out from the back).

The network science group (you can see me peeking out from the back, just left of center).

Anyway, it was a really fun and valuable experience, and the AMS pays for everything and a bag of chips (if by chips you mean more travel money to meet up with your collaborators and a ticket to the AMS Joint Mathematics Meeting the following January). I’m excited to blog about the work that come out of this, as network science is right up there with the coolest of topics in math and programming.

So if you’re eligible, keep an eye out for next year’s program.

Three Years Old, and an Idea for a Podcast

Happy birthday, Math ∩ Programming!

Today marks the end of the third year I’ve been writing Math ∩ Programming, and I’m excited to keep it going as I start my research career. In the last year I’ve started a secondary writing blog for some smaller, less technical bits, mostly to get thoughts out of my head. And while I could use this anniversary post to preview future Math ∩ Programming posts or review old favorites, I’ll instead share an idea that has been bouncing around my head for a few weeks. I’d love to hear your feedback in the comments.

I listen to podcasts and radio shows a lot, mostly storytelling and interviews. And they’re always bringing on these fancy-sounding people who write books on the New York Times Bestsellers list and who often have very interesting things to say. When discussing science they can often convey the ideas to the clueless listener, usually because it’s experimental science that’s naturally easy to understand (state the setup, state the results, hypothesize about the implications). But almost unilaterally there’s nothing substantive about math. All the mathematical content is popular math, how beautiful is \pi and such; math education, which I love to read and talk about but is common; or math history, which I’m not as interested in. And when there is some breakthrough, like Grigori Perelman solving a Millennium prize problem, the focus is entirely on the person and not the achievement. This isn’t specific to podcasts, but all news. I just happen to prefer my news in podcast form.

And so, aside from the myriad of excellent technical blogs by active researchers, what is there really that conveys the excitement I experience in theoretical computer science? There are publications like the ACM SIGACT monthly newsletter, which has a ton of book reviews and a handful of technical columns. Unfortunately it’s hidden behind a paywall, which basically immediately excludes it from being accessed by anyone not already embedded deep in academia. That being said it often has really interesting pieces like a poll by Bill Gasarch (2002, 2012) of researchers and their opinions on P vs NP. It’s really interesting to see just how much people differ on their desire to see other parts of mathematics incorporated into its resolution.

So if you don’t want to pay the ACM for a monthly newsletter, what can you do? Many of these ideas and opinions don’t exist in textbooks, and textbooks can be dry and bad at conveying why things are interesting or exciting. There are abstruse technical papers that you have to finish a graduate degree before you can even parse what’s being said. And then there are talks, which vary in quality almost as much as prose in technical papers do.

I recently came across a paper by Ryan Williams, a prominent researcher in circuit complexity. Roughly, when you study circuit complexity you try to understand which problems provably require big circuits to solve, and you study those proof techniques. It sounds boring but it’s interesting for three reasons: it’s extremely hard, there are many “embarrassing” open problems, and many of these problems imply wonderful things like P \neq NP. I actually get really excited by circuit complexity.

Anyway, this paper was titled “A Casual Tour Around a Circuit Complexity Bound,” in which Ryan reflects on the path which led him to one of the biggest breakthroughs in the last five years in circuit complexity. His writing is more or less informal (it was published in the SIGACT newsletter, though I had to access it through arXiv), and it focuses heavily on the big picture. It struck me as mostly how to think about circuit complexity. This kind of thing is truly invaluable for a graduate student and anyone, I imagine, trying to learn more about circuit complexity. Honestly, I’d love to see more of this in academic literature. Often papers are expanded from relatively simple principles into a mess of technical details, and reversing this process is slow and difficult.

But even besides these huge breakthroughs there are often really great ways to explain new problems and solutions. For example, this paper of Andrew Drucker, titled “High-Confidence Predictions under Adversarial Uncertainty,” starts with a really easy to understand setup:

A frog wants to cross the road at some fixed location, to get to a nice pond. But she is concerned about cars. It takes her a minute to cross the road, and if a car passes during that time, she will be squashed. However, this is no ordinary frog. She is extremely patient, and happy to wait any finite number of steps to cross the road. What’s more, she can observe and remember how many cars have passed, as well as when they passed. She can follow any algorithm to determine when to cross the road based on what she has seen so far, although her senses aren’t keen enough to detect a car before it arrives…

[Even if we assume the cars arrive according to a fixed probability distribution,] the frog may not have a detailed idea of how the cars are generated. It may be that the frog merely knows or conjectures some constraint obeyed by the car-stream. We then ask whether there exists a strategy which gets the frog safely across the road (at least, with sufficiently high probability), for any car-stream obeying the constraint.

This kind of story is better than coffee at keeping people awake during talks!

And so, I have been thinking a lot about what a podcast about theoretical computer science might entail. I imagine it going something like this: every episode is a half-hour conversation with a prominent researcher. The discussion would cover something about past work, something about future ideas of what’s important and a high level idea of the burgeoning techniques, and overarching questions about how one approaches research. Computer science is particularly interesting because most graduate students know enough to start working on open problems in their first year (so the topics are more accessible than, say, algebraic geometry), and because basically all of the theorems with names are named after people still active in the research community. Moreover the format of a podcast would require the interviewees to phrase their research in a way that doesn’t require a chalkboard or notation.

The hope is to popularize theoretical computer science by assuming some modest level of technical proficiency and to give access to anyone who wants to listen; say, an advanced undergraduate, an early graduate student, or a redditor who likes to argue about the Turing test. Moreover, if I were to actually run such a podcast, I could fill the listener in with additional details in a preface. The podcast could be a resource for undergraduates who want to explore the landscape of research topics before applying to graduate school, for graduate students who want to learn to think like a researcher or hear the variety of views out there, for researchers to advertise their favorite topics and get people to read/cite their papers, and for me to meet and interact with all of these great people. My advisor seems to know everyone and their lemmas, and more and more I’ve been finding myself wanting this as well (I’m just so terrible with names!). I’ve had enough conversations with researchers to know they have heaps of interesting things to say. And I travel enough to conferences and workshops. I imagine it wouldn’t be too hard to orchestrate a 30-minute conversation with one or two people. I’d just have to work up the courage to ask :)

What do you think of the idea? Would you listen to a theory podcast? Do you have a good idea for a catchy name? Are you a theory researcher who would like to have a conversation with me the next time we’re at a conference together? *fingers crossed* I’d love to hear from you.

Linear Programming and the Most Affordable Healthy Diet — Part 1

Optimization is by far one of the richest ways to apply computer science and mathematics to the real world. Everybody is looking to optimize something: companies want to maximize profits, factories want to maximize efficiency, investors want to minimize risk, the list just goes on and on. The mathematical tools for optimization are also some of the richest mathematical techniques. They form the cornerstone of an entire industry known as operations research, and advances in this field literally change the world.

The mathematical field is called combinatorial optimization, and the name comes from the goal of finding optimal solutions more efficiently than an exhaustive search through every possibility. This post will introduce the most central problem in all of combinatorial optimization, known as the linear program. Even better, we know how to efficiently solve linear programs, so in future posts we’ll write a program that computes the most affordable diet while meeting the recommended health standard.

Generalizing a Specific Linear Program

Most optimization problems have two parts: an objective function, the thing we want to maximize or minimize, and constraints, rules we must abide by to ensure we get a valid solution. As a simple example you may want to minimize the amount of time you spend doing your taxes (objective function), but you certainly can’t spend a negative amount of time on them (a constraint).

The following more complicated example is the centerpiece of this post. Most people want to minimize the amount of money spent on food. At the same time, one needs to maintain a certain level of nutrition. For males ages 19-30, the United States National Institute for Health recommends 3.7 liters of water per day, 1,000 milligrams of calcium per day, 90 milligrams of vitamin C per day, etc.

We can set up this nutrition problem mathematically, just using a few toy variables. Say we had the option to buy some combination of oranges, milk, and broccoli. Some rough estimates [1] give the following content/costs of these foods. For 0.272 USD you can get 100 grams of orange, containing a total of 53.2mg of calcium, 40mg of vitamin C, and 87g of water. For 0.100 USD you can get 100 grams of whole milk, containing 276mg of calcium, 0mg of vitamin C, and 87g of water. Finally, for 0.381 USD you can get 100 grams of broccoli containing 47mg of calcium, 89.2mg of vitamin C, and 91g of water. Here’s a table summarizing this information:

Nutritional content and prices for 100g of three foods

Food         calcium(mg)     vitamin C(mg)      water(g)   price(USD/100g)
Broccoli     47              89.2               91         0.381
Whole milk   276             0                  87         0.100
Oranges      40              53.2               87         0.272

Some observations: broccoli is more expensive but gets the most of all three nutrients, whole milk doesn’t have any vitamin C but gets a ton of calcium for really cheap, and oranges are a somewhere in between. So you could probably tinker with the quantities and figure out what the cheapest healthy diet is. The problem is what happens when we incorporate hundreds or thousands of food items and tens of nutrient recommendations. This simple example is just to help us build up a nice formality.

So let’s continue doing that. If we denote by b the number of 100g units of broccoli we decide to buy, and m the amount of milk and r the amount of oranges, then we can write the daily cost of food as

\displaystyle \text{cost}(b,m,r) = 0.381 b + 0.1 m + 0.272 r

In the interest of being compact (and again, building toward the general linear programming formulation) we can extract the price information into a single cost vector c = (0.381, 0.1, 0.272), and likewise write our variables as a vector x = (b,m,r). We’re implicitly fixing an ordering on the variables that is maintained throughout the problem, but the choice of ordering doesn’t matter. Now the cost function is just the inner product (dot product) of the cost vector and the variable vector \left \langle c,x \right \rangle. For some reason lots of people like to write this as c^Tx, where c^T denotes the transpose of a matrix, and we imagine that c and x are matrices of size 3 \times 1. I’ll stick to using the inner product bracket notation.

Now for each type of food we get a specific amount of each nutrient, and the sum of those nutrients needs to be bigger than the minimum recommendation. For example, we want at least 1,000 mg of calcium per day, so we require that 1000 \leq 47b + 276m + 40r. Likewise, we can write out a table of the constraints by looking at the columns of our table above.

\displaystyle \begin{matrix} 91b & + & 87m & + & 87r & \geq & 3700 & \text{(water)}\\ 47b & + & 276m & + & 40r & \geq & 1000 & \text{(calcium)} \\ 89.2b & + & 0m & + & 53.2r & \geq & 90 & \text{(vitamin C)} \end{matrix}

In the same way that we extracted the cost data into a vector to separate it from the variables, we can extract all of the nutrient data into a matrix A, and the recommended minimums into a vector v. Traditionally the letter b is used for the minimums vector, but for now we’re using b for broccoli.

A = \begin{pmatrix} 91 & 87 & 87 \\ 47 & 276 & 40 \\ 89.2 & 0 & 53.2 \end{pmatrix}

v = \begin{pmatrix} 3700 \\ 1000 \\ 90 \end{pmatrix}

And now the constraint is that Ax \geq v, where the \geq means “greater than or equal to in every coordinate.” So now we can write down the more general form of the problem for our specific matrices and vectors. That is, our problem is to minimize \left \langle c,x \right \rangle subject to the constraint that Ax \geq v. This is often written in offset form to contrast it with variations we’ll see in a bit:

\displaystyle \text{minimize} \left \langle c,x \right \rangle \\ \text{subject to the constraint } Ax \geq v

In general there’s no reason you can’t have a “negative” amount of one variable. In this problem you can’t buy negative broccoli, so we’ll add the constraints to ensure the variables are nonnegative. So our final form is

\displaystyle \text{minimize} \left \langle c,x \right \rangle \\ \text{subject to } Ax \geq v \\ \text{and } x \geq 0

In general, if you have an m \times n matrix A, a “minimums” vector v \in \mathbb{R}^m, and a cost vector c \in \mathbb{R}^n, the problem of finding the vector x that minimizes the cost function while meeting the constraints is called a linear programming problem or simply a linear program.

To satiate the reader’s burning curiosity, the solution for our calcium/vitamin C problem is roughly x = (1.01, 41.47, 0). That is, you should have about 100g of broccoli and 4.2kg of milk (like 4 liters), and skip the oranges entirely. The daily cost is about 4.53 USD. If this seems awkwardly large, it’s because there are cheaper ways to get water than milk.

100-grams-broccoli

100g of broccoli (image source: 100-grams.blogspot.com)

[1] Water content of fruits and veggiesFood costs in March 2014 in the midwest, and basic known facts about the water density/nutritional content of various foods.

Duality

Now that we’ve seen the general form a linear program and a cute example, we can ask the real meaty question: is there an efficient algorithm that solves arbitrary linear programs? Despite how widely applicable these problems seem, the answer is yes!

But before we can describe the algorithm we need to know more about linear programs. For example, say you have some vector x which satisfies your constraints. How can you tell if it’s optimal? Without such a test we’d have no way to know when to terminate our algorithm. Another problem is that we’ve phrased the problem in terms of minimization, but what about problems where we want to maximize things? Can we use the same algorithm that finds minima to find maxima as well?

Both of these problems are neatly answered by the theory of duality. In mathematics in general, the best way to understand what people mean by “duality” is that one mathematical object uniquely determines two different perspectives, each useful in its own way. And typically a duality theorem provides one with an efficient way to transform one perspective into the other, and relate the information you get from both perspectives. A theory of duality is considered beautiful because it gives you truly deep insight into the mathematical object you care about.

In linear programming duality is between maximization and minimization. In particular, every maximization problem has a unique “dual” minimization problem, and vice versa. The really interesting thing is that the variables you’re trying to optimize in one form correspond to the contraints in the other form! Here’s how one might discover such a beautiful correspondence. We’ll use a made up example with small numbers to make things easy.

So you have this optimization problem

\displaystyle \begin{matrix}  \text{minimize} & 4x_1+3x_2+9x_3 & \\  \text{subject to} & x_1+x_2+x_3 & \geq 6 \\  & 2x_1+x_3 & \geq 2 \\  & x_2+x_3 & \geq 1 & \\ & x_1,x_2,x_3 & \geq 0 \end{matrix}

Just for giggles let’s write out what A and c are.

\displaystyle A = \begin{pmatrix} 1 & 1 & 1 \\ 2 & 0 & 1 \\ 0 & 1 & 1 \end{pmatrix}, c = (4,3,9), v = (6,2,1)

Say you want to come up with a lower bound on the optimal solution to your problem. That is, you want to know that you can’t make 4x_1 + 3x_2 + 9x_3 smaller than some number m. The constraints can help us derive such lower bounds. In particular, every variable has to be nonnegative, so we know that 4x_1 + 3x_2 + 9x_3 \geq x_1 + x_2 + x_3 \geq 6, and so 6 is a lower bound on our optimum. Likewise,

\displaystyle \begin{aligned}4x_1+3x_2+9x_3 & \geq 4x_1+4x_3+3x_2+3x_3 \\ &=2(2x_1 + x_3)+3(x_2+x_3) \\ & \geq 2 \cdot 2 + 3 \cdot 1 \\ &=7\end{aligned}

and that’s an even better lower bound than 6. We could try to write this approach down in general: find some numbers y_1, y_2, y_3 that we’ll use for each constraint to form

\displaystyle y_1(\text{constraint 1}) + y_2(\text{constraint 2}) + y_3(\text{constraint 3})

To make it a valid lower bound we need to ensure that the coefficients of each of the x_i are smaller than the coefficients in the objective function (i.e. that the coefficient of x_1 ends up less than 4). And to make it the best lower bound possible we want to maximize what the right-hand-size of the inequality would be: y_1 6 + y_2 2 + y_3 1. If you write out these equations and the constraints you get our “lower bound” problem written as

\displaystyle \begin{matrix} \text{maximize} & 6y_1 + 2y_2 + y_3 & \\ \text{subject to} & y_1 + 2y_2 & \leq 4 \\ & y_1 + y_3 & \leq 3 \\ & y_1+y_2 + y_3 & \leq 9 \\ & y_1,y_2,y_3 & \geq 0 \end{matrix}

And wouldn’t you know, the matrix providing the constraints is A^T, and the vectors c and v switched places.

\displaystyle A^T = \begin{pmatrix} 1 & 2 & 0 \\ 1 & 0 & 1 \\ 1 & 1 & 1 \end{pmatrix}

This is no coincidence. All linear programs can be transformed in this way, and it would be a useful exercise for the reader to turn the above maximization problem back into a minimization problem by the same technique (computing linear combinations of the constraints to make upper bounds). You’ll be surprised to find that you get back to the original minimization problem! This is part of what makes it “duality,” because the dual of the dual is the original thing again. Often, when we fix the “original” problem, we call it the primal form to distinguish it from the dual form. Usually the primal problem is the one that is easy to interpret.

(Note: because we’re done with broccoli for now, we’re going to use b to denote the constraint vector that used to be v.)

Now say you’re given the data of a linear program for minimization, that is the vectors c, b and matrix A for the problem, “minimize \left \langle c, x \right \rangle subject to Ax \geq b; x \geq 0.” We can make a general definition: the dual linear program is the maximization problem “maximize \left \langle b, y \right \rangle subject to A^T y \leq c, y \geq 0.” Here y is the new set of variables and the superscript T denotes the transpose of the matrix. The constraint for the dual is often written y^T A \leq c^T, again identifying vectors with a single-column matrices, but I find the swamp of transposes pointless and annoying (why do things need to be columns?).

Now we can actually prove that the objective function for the dual provides a bound on the objective function for the original problem. It’s obvious from the work we’ve done, which is why it’s called the weak duality theorem.

Weak Duality Theorem: Let c, A, b be the data of a linear program in the primal form (the minimization problem) whose objective function is \left \langle c, x \right \rangle. Recall that the objective function of the dual (maximization) problem is \left \langle b, y \right \rangle. If x,y are feasible solutions (satisfy the constraints of their respective problems), then

\left \langle b, y \right \rangle \leq \left \langle c, x \right \rangle

In other words, the maximum of the dual is a lower bound on the minimum of the primal problem and vice versa. Moreover, any feasible solution for one provides a bound on the other.

Proof. The proof is pleasingly simple. Just inspect the quantity \left \langle A^T y, x \right \rangle = \left \langle y, Ax \right \rangle. The constraints from the definitions of the primal and dual give us that

\left \langle y, b \right \rangle \leq \left \langle y, Ax \right \rangle = \left \langle A^Ty, x \right \rangle \leq \left \langle c,x \right \rangle

The inequalities follow from the linear algebra fact that if the u in \left \langle u,v \right \rangle is nonnegative, then you can only increase the size of the product by increasing the components of v. This is why we need the nonnegativity constraints.

In fact, the world is much more pleasing. There is a theorem that says the two optimums are equal!

Strong Duality Theorem: If there are any solutions x,y to the primal (minimization) problem and the dual (maximization) problem, respectively, then the two problems also have optimal solutions x^*, y^*, and two candidate solutions x^*, y^* are optimal if and only if they produce equal objective values \left \langle c, x^* \right \rangle = \left \langle y^*, b \right \rangle.

The proof of this theorem is a bit more convoluted than the weak duality theorem, and the key technique is a lemma of Farkas and its variations. See the second half of these notes for a full proof. The nice thing is that this theorem gives us a way to tell if an algorithm to solve linear programs is done: maintain a pair of feasible solutions to the primal and dual problems, improve them by some rule, and stop when the two solutions give equal objective values. The hard part, then, is finding a principled and guaranteed way to improve a given pair of solutions.

On the other hand, you can also prove the strong duality theorem by inventing an algorithm that provably terminates. We’ll see such an algorithm, known as the simplex algorithm in the next post. Sneak peek: it’s a lot like Gaussian elimination. Then we’ll use the algorithm (or an equivalent industry-strength version) to solve a much bigger nutrition problem.

In fact, you can do a bit better than the strong duality theorem, in terms of coming up with a stopping condition for a linear programming algorithm. You can observe that an optimal solution implies further constraints on the relationship between the primal and the dual problems. In particular, this is called the complementary slackness conditions, and they essentially say that if an optimal solution to the primal has a positive variable then the corresponding constraint in the dual problem must be tight (is an equality) to get an optimal solution to the dual. The contrapositive says that if some constraint is slack, or a strict inequality, then either the corresponding variable is zero or else the solution is not optimal. More formally,

Theorem (Complementary Slackness Conditions): Let A, c, b be the data of the primal form of a linear program, “minimize \left \langle c, x \right \rangle subject to Ax \geq b, x \geq 0.” Then x^*, y^* are optimal solutions to the primal and dual problems if any only if all of the following conditions hold.

  • x^*, y^* are both feasible for their respective problems.
  • Whenever x^*_i > 0 the corresponding constraint A^T_i y^* = c_i is an equality.
  • Whenever y^*_j > 0 the corresponding constraint A_j x^* = b_j is an equality.

Here we denote by M_i the i-th row of the matrix M and v_i to denote the i-th entry of a vector. Another way to write the condition using vectors instead of English is

\left \langle x^*, A^T y^* - c \right \rangle = 0
\left \langle y^*, Ax^* - b \right \rangle

The proof follows from the duality theorems, and just involves pushing around some vector algebra. See section 6.2 of these notes.

One can interpret complementary slackness in linear programs in a lot of different ways. For us, it will simply be a termination condition for an algorithm: one can efficiently check all of these conditions for the nonzero variables and stop if they’re all satisfied or if we find a variable that violates a slackness condition. Indeed, in more mature optimization analyses, the slackness condition that is more egregiously violated can provide evidence for where a candidate solution can best be improved. For a more intricate and detailed story about how to interpret the complementary slackness conditions, see Section 4 of these notes by Joel Sobel.

Finally, before we close we should note there are geometric ways to think about linear programming. I have my preferred visualization in my head, but I have yet to find a suitable animation on the web that replicates it. Here’s one example in two dimensions. The set of constraints define a convex geometric region in the plane

The constraints define a convex area of "feasible solutions." Image source: Wikipedia.

The constraints define a convex area of “feasible solutions.” Image source: Wikipedia.

Now the optimization function f(x) = \left \langle c,x \right \rangle is also a linear function, and if you fix some output value y = f(x) this defines a line in the plane. As y changes, the line moves along its normal vector (that is, all these fixed lines are parallel). Now to geometrically optimize the target function, we can imagine starting with the line f(x) = 0, and sliding it along its normal vector in the direction that keeps it in the feasible region. We can keep sliding it in this direction, and the maximum of the function is just the last instant that this line intersects the feasible region. If none of the constraints are parallel to the family of lines defined by f, then this is guaranteed to occur at a vertex of the feasible region. Otherwise, there will be a family of optima lying anywhere on the line segment of last intersection.

In higher dimensions, the only change is that the lines become affine subspaces of dimension n-1. That means in three dimensions you’re sliding planes, in four dimensions you’re sliding 3-dimensional hyperplanes, etc. The facts about the last intersection being a vertex or a “line segment” still hold. So as we’ll see next time, successful algorithms for linear programming in practice take advantage of this observation by efficiently traversing the vertices of this convex region. We’ll see this in much more detail in the next post.

Until then!