|Making Hybrid Images||Neural Networks
|Bezier Curves and Picasso||Computing Homology||Probably Approximately
Correct – A Formal Theory
Problem: Alice chooses a secret polynomial with nonnegative integer coefficients. Bob wants to discover this polynomial by querying Alice for the value of for some integer of Bob’s choice. What is the minimal number of queries Bob needs to determine exactly?
Solution: Two queries. The first is , and if we call , then the second query is .
To someone who is familiar with polynomials, this may seem shocking, and I’ll explain why it works in a second. After all, it’s very easy to prove that if Bob gives Alice all of his queries at the same time (if the queries are not adaptive), then it’s impossible to discover what is using fewer than queries. This is due to a fact called polynomial interpolation, which we’ve seen on this blog before in the context of secret sharing. Specifically, there is a unique single-variable degree polynomial passing through points (with distinct -values). So if you knew the degree of , you could determine it easily. But Bob doesn’t know the degree of the polynomial, and there’s no way he can figure it out without adaptive queries! Indeed, if Bob tries and gives a set of queries, Alice could have easily picked a polynomial of degree . So it’s literally impossible to solve this problem without adaptive queries.
The lovely fact is that once you allow adaptiveness, the number of queries you need doesn’t even depend on the degree of the secret polynomial!
Okay let’s get to the solution. It was crucial that our polynomial had nonnegative integer coefficients, because we’re going to do a tiny bit of number theory. Let . First, note that is exactly the sum of the coefficients , and in particular is larger than any single coefficient. So call this , and query . This gives us a number of the form
And because is so big, we can compute easily by computing . Now set , and this has the form . We can compute modulus again to get , and repeat until we have all the coefficients. We’ll stop once we get a that is zero.
As a small technical note, this is a polynomial-time algorithm in the number of bits needed to write down . So this demonstrates the power of adaptive queries: we get from something which is uncomputable with any number of queries to something which is efficiently computable with a constant number of queries.
The obvious follow-up question is: can you come up with an efficient algorithm if we allow the coefficients to be negative integers?
One of the most interesting questions posed in the last thirty years of computer science is to ask how much “information” must be communicated between two parties in order for them to jointly compute something. One can imagine these two parties living on distant planets, so that the cost of communicating any amount of information is very expensive, but each person has an integral component of the answer that the other does not.
Since this question was originally posed by Andrew Yao in 1979, it has led to a flurry of applications in many areas of mathematics and computer science. In particular it has become a standard tool for proving lower bounds in many settings such as circuit design and streaming algorithms. And if there’s anything theory folks love more than a problem that can be solved by an efficient algorithm, it’s a proof that a problem cannot be solved by any efficient algorithm (that’s what I mean by “lower bound”).
Despite its huge applicability, the basic results in this area are elementary. In this post we’ll cover those basics, but once you get past these basic ideas and their natural extensions you quickly approach the state of the art and open research problems. Attempts to tackle these problems in recent years have used sophisticated techniques in Fourier analysis, Ramsey theory, and geometry. This makes it a very fun and exciting field.
As a quick side note before we start, the question we’re asking is different from the one of determining the information content of a specific message. That is the domain of information theory, which was posed (and answered) decades earlier. Here we’re trying to determine the complexity of a problem, where more complex messages require more information about their inputs.
The Basic Two-Player Model
The most basic protocol is simple enough to describe over a dinner table. Alice and Bob each have one piece of information , respectively, say they each have a number. And together they want to compute some operation that depends on both their inputs, for example whether . But in the beginning Alice has access only to her number , and knows nothing about . So Alice sends Bob a few bits. Depending on the message Bob computes something and replies, and this repeats until they have computed an answer. The question is: what is the minimum number of bits they need to exchange in order for both of them to be able to compute the right answer?
There are a few things to clarify here: we’re assuming that Alice and Bob have agreed on a protocol for sending information before they ever saw their individual numbers. So imagine ten years earlier Alice and Bob were on the same planet, and they agreed on the rules they’d follow for sending/replying information once they got their numbers. In other words, we’re making a worst-case assumption on Alice and Bob’s inputs, and as usual it will be measured as a function of , the lengths of their inputs. Then we take a minimum (asymptotically) over all possible protocols they could follow, and this value is the “communication complexity” of the problem. Computing the exact communication complexity of a given problem is no simple task, since there’s always the nagging question of whether there’s some cleverer protocol than the one you came up with. So most of the results are bounds on the communication complexity of a problem.
Indeed, we can give our first simple bound for the “ greater than ” problem we posed above. Say the strings both have bits. What Alice does is send her entire string to Bob, and Bob then computes the answer and sends the answer bit back to Alice. This requires bits of communication. This proves that the communication complexity of “” is bounded from above by . A much harder question is, can we do any better?
To make any progress on upper or lower bounds we need to be a bit more formal about the communication model. Basically, the useful analysis happens when the players alternate sending single bits, and this is only off by small constant factors from a more general model. This is the asymptotic analysis, that we only distinguish between things like linear complexity versus sublinear options like or or even constant complexity . Indeed, the protocol we described for is the stupidest possible protocol for the problem, and it’s actually valid for any problem. For this problem it happens to be optimal, but we’re just trying to emphasize that nontrivial bounds are all sub-linear in the size of the inputs.
On to the formal model.
Definition: A player is a computationally unbounded Turing machine.
And we really mean unbounded. Our players have no time or space constraints, and if they want they can solve undecidable problems like the halting problem or computing Kolmogorov complexity. This is to emphasize that the critical resource is the amount of communication between players. Moreover, it gives us a hint that lower bounds in this model won’t come form computational intractability, but instead will be “information-theoretic.”
Definition: Let be the set of all binary strings. A communication protocol is a pair of functions .
The input to these functions should be thought of as follows: is the player’s secret input and is the communication history so far. The output is the single bit that they will send in that round (which can be determined by the length of since only one bit is sent in each round). The protocol then runs by having Alice send to Bob, then Bob replies with , Alice continues with , and so on. We implicitly understand that the content of a communication protocol includes a termination condition, but we’ll omit this from the notation. We call the length of the protocol the number of rounds.
Definition: A communication protocol is said to be valid for a boolean function if for all strings , the protocol for terminates on some round with if and only if .
So to define the communication complexity, we let the function be the maximum length of the protocol when run on strings of length (the worst-case for a given input size). Then the communication complexity of a function is the minimum of over all valid protocols . In symbols,
We will often abuse the notation by writing the communication complexity of a function as , understanding that it’s measured asymptotically as a function of .
Matrices and Lower Bounds
Let’s prove a lower bound, that to compute the equality function you need to send a linear number of bits in the worst case. In doing this we’ll develop a general algebraic tool.
So let’s write out the function as a binary matrix in the following way. Write all inputs of length in some fixed order along the rows and columns of the matrix, and let entry be the value of . For example, the 6-bit function which computes whether the majority of the two player’s bits are ones looks like this:
The key insight to remember is that if the matrix of a function has a nice structure, then one needs very little communication to compute it. Let’s see why.
Say in the first round the row player sends a bit . This splits the matrix into two submatrices by picking the rows of to be those inputs for which the row player sends a , and likewise for with . If you’re willing to rearrange the rows of the matrix so that and stack on top of each other, then this splits the matrix into two rectangles. Now we can switch to the column player and see which bit he sends in reply to each of the possible choices for (say he sends back ). This separately splits each of into two subrectangles corresponding to which inputs for the column player make him send the specific value of . Continuing in this fashion we recurse until we find a submatrix consisting entirely of ones or entirely of zeros, and then we can say with certainty what the value of the function is.
It’s difficult to visualize because every time we subdivide we move around the rows and columns within the submatrix corresponding to the inputs for each player. So the following would be a possible subdivision of an 8×8 matrix (with the values in the rectangles denoting which communicated bits got you there), but it’s sort of a strange one because we didn’t move the inputs around at all. It’s just a visual aid.
If we do this for steps we get subrectangles. A crucial fact is that any valid communication protocol for a function has to give a subdivision of the matrix where all the rectangles are constant. or else there would be two pairs of inputs , which are labeled identically by the communication protocol, but which have different values under .
So naturally one expects the communication complexity of would require as many steps as there are steps in the best decomposition, that is, the decomposition with the fewest levels of subdivision. Indeed, we’ll prove this and introduce some notation to make the discourse less clumsy.
Definition: For an matrix , a rectangle is a submatrix where . A rectangle is called monochromatic if all entires in the corresponding submatrix are the same. A monochromatic tiling of is a partition of into disjoint monochromatic rectangles. Define to be the minimum number of rectangles in any monochromatic tiling of .
As we said, if there are steps in a valid communication protocol for , then there are rectangles in the corresponding monochromatic tiling of . Here is an easy consequence of this.
Proposition: If has communication complexity , then there is a monochromatic tiling of with at most rectangles. In particular, .
Proof. Pick any protocol that achieves the communication complexity of , and apply the process we described above to subdivide . This will take exactly , and produce no more than rectangles.
This already gives us a bunch of theorems. Take the EQ function, for example. Its matrix is the identity matrix, and it’s not hard to see that every monochromatic tiling requires rectangles, one for each entry of the diagonal. I.e., . But we already know that one player can just send all his bits, so actually . Now it’s not always so easy to compute . The impressive thing to do is to use efficiently computable information about to give bounds on and hence on . So can we come up with a better lower bound that depends on something we can compute? The answer is yes.
Theorem: For every function , .
Proof. This just takes some basic linear algebra. One way to think of the rank of a matrix is as the smallest way to write as a linear combination of rank 1 matrices (smallest as in, the smallest number of terms needed to do this). The theorem is true no matter which field you use to compute the rank, although in this proof and in the rest of this post we’ll use the real numbers.
If you give me a monochromatic tiling by rectangles, I can view each rectangle as a matrix whose rank is at most one. If the entries are all zeros then the rank is zero, and if the entries are all ones then (using zero elsewhere) this is by itself a rank 1 matrix. So adding up these rectangles as separate components gives me an upper bound on the rank of . So the minimum way to do this is also an upper bound on the rank of .
Now computing something like is even easier, because the rank of is just .
There are other techniques to show lower bounds that are stronger than the rank and tiling method (because they imply the rank and tiling method). See this survey for a ton of details. But I want to discuss upper bounds a bit, because the central open conjecture in communication complexity is an upper bound.
The Log-Rank Conjecture: There is a universal constant , such that for all , the communication complexity .
All known examples satisfy the conjecture, but unfortunately the farthest progress toward the conjecture is still exponentially worse than the conjecture’s statement. In 1997 the record was due to Andrei Kotlov who proved that . It was not until 2013 that any (unconditional) improvements were made to this, when Shachar Lovett proved that .
The interested reader can check out this survey of Shachar Lovett from earlier this year (2014) for detailed proofs of these theorems and a discussion of the methods. I will just discuss one idea from this area that ties in nicely with our discussion: which is that finding an efficient communication protocol for a low-rank function reduces to finding a large monochromatic rectangle in .
Theorem [Nisan-Wigderson 94]: Let be a function. Suppose that for any function , we can find a monochromatic rectangle of size where . Then any such is computable by a deterministic protocol with communication complexity.
Just to be concrete, this says that if is polylogarithmic, then finding these big rectangles implies a protocol also with polylogarithmic complexity. Since the complexity of the protocol is a function of alone, the log-rank conjecture follows as a consequence. The best known results use the theorem for larger for some , which gives communication complexity also .
The proof of the theorem is detailed, but mostly what you’d expect. You take your function, split it up into the big monochromatic rectangle and the other three parts. Then you argue that when you recurse to one of the other three parts, either the rank is cut in half, or the size of the matrix is much smaller. In either case, you can apply the theorem once again. Then you bound the number of leaves in the resulting protocol tree by looking at each level where the rank has dropped to . For the full details, see page 4 of the Shachar survey.
Multiple Players and More
In the future we’ll cover some applications of communication complexity, many of which are related to computing in restricted models such as parallel computation and streaming computation. For example, in parallel computing you often have processors which get arbitrary chunks of data as input and need to jointly compute something. Lower bounds on the communication complexity can help you prove they require a certain amount of communication in order to do that.
But in these models there are many players. And the type of communication matters: it can be point-to-point or broadcast, or something more exotic like MapReduce. So before we can get to these applications we need to define and study the appropriate generalizations of communication complexity to multiple interacting parties.
I recently wrapped up a fun paper with my coauthors Ben Fish, Adam Lelkes, Lev Reyzin, and Gyorgy Turan in which we analyzed the computational complexity of a model of the popular MapReduce framework. Check out the preprint on the arXiv.
As usual I’ll give a less formal discussion of the research here, and because the paper is a bit more technically involved than my previous work I’ll be omitting some of the more pedantic details. Our project started after Ben Moseley gave an excellent talk at UI Chicago. He presented a theoretical model of MapReduce introduced by Howard Karloff et al. in 2010, and discussed his own results on solving graph problems in this model, such as graph connectivity. You can read Karloff’s original paper here, but we’ll outline his model below.
Basically, the vast majority of the work on MapReduce has been algorithmic. What I mean by that is researchers have been finding more and cleverer algorithms to solve problems in MapReduce. They have covered a huge amount of work, implementing machine learning algorithms, algorithms for graph problems, and many others. In Moseley’s talk, he posed a question that caught our eye:
Is there a constant-round MapReduce algorithm which determines whether a graph is connected?
After we describe the model below it’ll be clear what we mean by “solve” and what we mean by “constant-round,” but the conjecture is that this is impossible, particularly for the case of sparse graphs. We know we can solve it in a logarithmic number of rounds, but anything better is open.
In any case, we started thinking about this problem and didn’t make much progress. To the best of my knowledge it’s still wide open. But along the way we got into a whole nest of more general questions about the power of MapReduce. Specifically, Karloff proved a theorem relating MapReduce to a very particular class of circuits. What I mean is he proved a theorem that says “anything that can be solved in MapReduce with so many rounds and so much space can be solved by circuits that are yae big and yae complicated, and vice versa.
But this question is so specific! We wanted to know: is MapReduce as powerful as polynomial time, our classical notion of efficiency (does it equal P)? Can it capture all computations requiring logarithmic space (does it contain L)? MapReduce seems to be somewhere in between, but it’s exact relationship to these classes is unknown. And as we’ll see in a moment the theoretical model uses a novel communication model, and processors that never get to see the entire input. So this led us to a host of natural complexity questions:
- What computations are possible in a model of parallel computation where no processor has enough space to store even one thousandth of the input?
- What computations are possible in a model of parallel computation where processor’s can’t request or send specific information from/to other processors?
- How the hell do you prove that something can’t be done under constraints of this kind?
- How do you measure the increase of power provided by giving MapReduce additional rounds or additional time?
These questions are in the domain of complexity theory, and so it makes sense to try to apply the standard tools of complexity theory to answer them. Our paper does this, laying some brick for future efforts to study MapReduce from a complexity perspective.
In particular, one of our theorems is the following:
Theorem: Any problem requiring sublogarithmic space can be solved in MapReduce in two rounds.
This theorem is nice for two reasons. First is it’s constructive. If you give me a problem and I know it classically takes less than logarithmic space, then this gives an automatic algorithm to implement it in MapReduce, and if you were so inclined you could even automate the process of translating a classical algorithm to a MapReduce algorithm (it’s not pretty, but it’s possible).
Our other results are a bit more esoteric. We relate the questions about MapReduce to old, deep questions about computations that have simultaneous space and time bounds. In the end we give a (conditional, nonconstructive) answer to question 4 above, which I’ll sketch without getting too deep in the details.
So let’s start with the model of Karloff et al., which they named “MRC” for MapReduce Class.
The MRC Model of Karloff et al.
MapReduce is a programming paradigm which is intended to make algorithm design for distributed computing easier. Specifically, when you’re writing massively distributed algorithms by hand, you spend a lot of time and effort dealing with really low-level questions. Like, what do I do if a processor craps out in the middle of my computation? Or, what’s the most efficient way to broadcast a message to all the processors, considering the specific topology of my network layout means the message will have to be forwarded? Note these are questions that have nothing to do with the actual problem you’re trying to solve.
So the designers of MapReduce took a step back, analyzed the class of problems they were most often trying to solve (turns out, searching, sorting, and counting), and designed their framework to handle all of the low-level stuff automatically. The catch is that the algorithm designer now has to fit their program into the MapReduce paradigm, which might be hard.
In designing a MapReduce algorithm, one has to implement two functions: a mapper and a reducer. The mapper takes a list of key-value pairs, and applies some operation to each element independently (just like the map function in most functional programming languages). The reducer takes a single key and a list of values for that key, and outputs a new list of values with the same key. And then the MapReduce protocol iteratively applies mappers and reducers in rounds to the input data. There are a lot of pictures of this protocol on the internet. Here’s one
An interesting part of this protocol is that the reducers never get to communicate with each other, except indirectly through the mappers in the next round. MapReduce handles the implied grouping and message passing, as well as engineering issues like fault-tolerance and load balancing. In this sense, the mappers and reducers are ignorant cogs in the machine, so it’s interesting to see how powerful MapReduce is.
The model of Karloff, Suri, and Vassilvitskii is a relatively straightforward translation of this protocol to mathematics. They make a few minor qualifications, though, the most important of which is that they encode a bound on the total space used. In particular, if the input has size , they impose that there is some for which every reducer uses space . Moreover, the number of reducers in any round is also bounded by .
We can formalize all of this with a few easy definitions.
Definition: An input to a MapReduce algorithm is a list of key-value pairs of total size .
For binary languages (e.g., you give me a binary string and you want to know if there are an odd number of 1’s), we encode a string as the list . This won’t affect any of our space bounds, because the total blowup from to is logarithmic.
Definition: A mapper is a Turing machine which accepts as input a single key-value pair and produces as output a list of key-value pairs .
Definition: A reducer is a Turing machine which accepts as input a key and a list of values and produces as output the same key and a new list of values .
Now we can describe the MapReduce protocol, i.e. what a program is and how to run it. I copied this from our paper because the notation is the same so far.
All we’ve done here is just describe the MapReduce protocol in mathematics. It’s messier than it is complicated. The last thing we need is the space bounds.
Definition: A language is in if there is a constant and a sequence of mappers and reducers such that for all the following is satisfied:
- Letting and , accepts if and only if .
- For all , use space and time.
- Each outputs keys in round .
To be clear, the first argument to is the round bound, and the second argument is the time bound. The last point implies the bound on the number of processors. Since we are often interested in a logarithmic number of rounds, we can define
So we can restate the question from the beginning of the post as,
Is graph connectivity in ?
Here there’s a bit of an issue with representation. You have to show that if you’re just given a binary string representing a graph, that you can translate that into a reasonable key-value description of a graph in a constant number of rounds. This is possible, and gives evidence that the key-value representation is without loss of generality for all intents and purposes.
A Pedantic Aside
If our goal is to compare MRC with classes like polynomial time (P) and logarithmic space (L), then the definition above has a problem. Specifically the definition above allows one to have an infinite list of reducers, where each one is potentially different, and the actual machine that is used depends on the input size. In formal terminology, MRC as defined above is a nonuniform model of computation.
Indeed, we give a simple proof that this is a problem by showing any unary language is in (which is where many of the MRC algorithms in recent years have been). For those who don’t know, this is a problem because you can encode versions of the Halting problem as a unary language, and the Halting problem is undecidable.
While this level of pedantry might induce some eye-rolling, it is necessary in order to make statements like “MRC is contained in P.” It’s simply not true for the nonuniform model above. To fix this problem we refined the MRC model to a uniform version. The details are not trivial, but also not that interesting. Check out the paper if you want to know exactly how we do it. It doesn’t have that much of an effect on the resulting analysis. The only important detail is that now we’re representing the entire MRC machine as a single Turing machine. So now, unlike before, any MRC machine can be written down as a finite string, and there are infinitely many strings representing a given MRC machine. We call our model MRC, and Karloff’s model nonuniform MRC.
You can also make randomized versions of MRC, but we’ll stick to the deterministic version here.
Now we can sketch the proof that sublogarithmic space is in . In fact, the proof is simpler for regular languages (equivalent to constant space algorithms) being in . The idea is as follows:
A regular language is one that can be decided by a DFA, say (a graph representing state transitions as you read a stream of bits). And the DFA is independent of the input size, so every mapper and reducer can have it encoded into them. Now what we’ll do is split the input string up in contiguous chunks of size among the processors (mappers can do this just fine). Now comes the trick. We have each reducer compute, for each possible state in , what the ending state would be if the DFA had started in state after processing their chunk of the input. So the output of reducer would be an encoding of a table of the form:
And the result would be a lookup table of intermediate computation results. So each reducer outputs their part of the table (which has constant size). Since there are only reducers, all of the table pieces can fit on a single reducer in the second round, and this reducer can just chain the functions together from the starting state and output the answer.
The proof for sublogarithmic space has the same structure, but is a bit more complicated because one has to account for the tape of a Turing machine which has size . In this case, you split up the tape of the Turing machine among the processors as well. And then you have to keep track of a lot more information. In particular, each entry of your function has to look like
“if my current state is A and my tape contents are B and the tape head starts on side C of my chunk of the input”
“the next time the tape head would leave my chunk, it will do so on side C’ and my state will be A’ and my tape contents will be B’.”
As complicated as it sounds, the size of one of these tables for one reducer is still less than for some . And so we can do the same trick of chaining the functions together to get the final answer. Note that this time the chaining will be adaptive, because whether the tape head leaves the left or right side will determine which part of the table you look at next. And moreover, we know the chaining process will stop in polynomial time because we can always pick a Turing machine to begin with that will halt in polynomial time (i.e., we know that L is contained in P).
The size calculations are also just large enough to stop us from doing the same trick with logarithmic space. So that gives the obvious question: is ? The second part of our paper shows that even weaker containments are probably very hard to prove, and they relate questions about MRC to the problem of L vs P.
This part of the paper is where we go much deeper into complexity theory than I imagine most of my readers are comfortable with. Our main theorem looks like this:
Theorem: Assume some conjecture is true. Then for every , there is are bigger such that the following hold:
This is a kind of hierarchy theorem that one often likes to prove in complexity theory. It says that if you give MRC enough extra rounds or time, then you’ll get strictly more computational power.
The conjecture we depend on is called the exponential time hypothesis (ETH), and it says that the 3-SAT problem can’t be solved in time for any . Actually, our theorem depends on a weaker conjecture (implied by ETH), but it’s easier to understand our theorems in the context of the ETH. Because 3-SAT is this interesting problem: we believe it’s time-intensive, and yet it’s space efficient because we can solve it in linear space given time. This time/space tradeoff is one of the oldest questions in all of computer science, but we still don’t have a lot of answers.
For instance, define to the the class of languages that can be decided by Turing machines that have simultaneous bounds of time and space. For example, we just said that , and there is a famous theorem that says that SAT is not in ; apparently this is the best we know. So clearly there are some very basic unsolved problems about TISP. An important one that we address in our paper is whether there are hierarchies in TISP for a fixed amount of space. This is the key ingredient in proving our hierarchy theorem for MRC. In particular here is an open problem:
Conjecture: For every two integers , the classes and are different.
We know this is true of time and space separately, i.e., that and . but for TISP we only know that you get more power if you increase both parameters.
So we prove a theorem that address this, but falls short in two respects: it depends on a conjecture like ETH, and it’s not for every .
Theorem: Suppose ETH holds, then for every there is a for which .
In words, this suggests that there are problems that need both time and space, but can be solved with linear space if you allow enough extra time. Since , this gives us the hierarchy we wanted.
To prove this we take 3-SAT and give it exponential padding so that it becomes easy enough to do in polynomial TISP (and it still takes linear space, in fact sublinear), but not so easy that you can do it in time. It takes some more work to get from this TISP hierarchy to our MRC hierarchy, but the details are a bit much for this blog. One thing I’d like to point out is that we prove that statements that are just about MRC have implications beyond MapReduce. In particular, the last corollary of our paper is the following.
Corollary: If , then L is different from P.
In other words, if you’re afforded a polynomial number of rounds in MRC, then showing that constant time per round is weaker than polynomial time is equivalently hard to separating L from P. The theorem is true because, as it turns out, , by simulating one step of a TM across polynomially many rounds. The proof is actually not that complicated (and doesn’t depend on ETH at all), and it’s a nice demonstration that studying MRC can have implications beyond parallel computing.
The other side of the coin is also interesting. Our first theorem implies the natural question of whether . We’d like to say that this would imply the separation of L from P, but we don’t quite get that. In particular, we know that
But at the same time we could live in a world where
It seems highly unlikely, but to the best of our knowledge none of our techniques prove this is not the case. If we could rule this out, then we could say that implies the separation of L and P. And note this would not depend on any conjectures.
Our paper has a list of open problems at the end. My favorite is essentially: how do we prove better lower bounds in MRC? In particular, it would be great if we could show that some problems need a lot of rounds simply because the communication model is too restrictive, and nobody has true random access to the entire input. For example, this is why we think graph connectivity needs a logarithmic number of rounds. But as of now nobody really knows how to prove it, and it seems like we’ll need some new and novel techniques in order to do it. I only have the wisps of ideas in that regard, and it will be fun to see which ones pan out.
Until next time!
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.”
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.
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.
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 with specific properties I won’t detail here like “smoothness” and “finiteness.” And the building blocks are the complex exponential functions
where 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 , 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 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 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:
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 and center is the function
It looks like this
In particular, at zero the function is 1 and it gradually drops to zero as you get farther away. The parameter controls the rate at which it vanishes, and in the picture above the center is set to .
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 , and if we’re doing a “high-pass filter” we multiply by .
Before we get to the code, here’s an example of a low-pass filter. First, take this image of Marilyn Monroe
Now compute its Fourier transform
Apply the low-pass filter
And reverse the Fourier transform to get an image
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
And if we compute the average of these two images, we basically get back to the original.
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).
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.
You can see some of the other hybrid images Oliva et al constructed over at their web gallery.
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!