Featured Posts

Making Hybrid Images Neural Networks
and Backpropagation
Elliptic Curves
and Cryptography
bezier-dog-image-small circle-wedge-sphere baby_programmers
Bezier Curves and Picasso Computing Homology Probably Approximately
Correct – A Formal Theory
of Learning

For mathematicians, = does not mean equality

Every now and then I hear some ridiculous things about the equals symbol. Some large subset of programmers—perhaps related to functional programmers, perhaps not—seem to think that = should only and ever mean “equality in the mathematical sense.” The argument usually goes,

Functional programming gives us back that inalienable right to analyze things by using mathematics. Never again need we bear the burden of that foul mutant x = x+1! No novice programmer—nay, not even a mathematician!—could comprehend such flabbergastery. Tis a pinnacle of confusion!

It’s ironic that so much of the merits or detriment of the use of = is based on a veiled appeal to the purity of mathematics. Just as often software engineers turn the tables, and any similarity to mathematics is decried as elitist jibber jabber (Such an archaic and abstruse use of symbols! Oh no, big-O!).

In fact, equality is more rigorously defined in a programming language than it will ever be in mathematics. Even in the hottest pits of software hell, where there’s = and == and ===, throwing in ==== just to rub salt in the wound, each operator gets its own coherent definition and documentation. Learn it once and you’ll never go astray.

Not so in mathematics—oh yes, hide your children from the terrors that lurk. In mathematics equality is little more than a stand-in for the word “is,” oftentimes entirely dependent on context. Now gather round and listen to the tale of the true identities of the masquerader known as =.

Let’s start with some low-hanging fruit, the superficial concerns.

\displaystyle \sum_{i=1}^n i^2 + 3

If = were interpreted literally, i would be “equal” to 1, and “equal” to 2, and I’d facetiously demand 1 = 2. Aha! Where is your Gauss now?! But seriously, this bit of notation shows that mathematics has both expressions with scope and variables that change their value over time. And the \sum use for notation was established by Euler, long before algorithms jumped from logic to computers to billionaire Senate testimonies.

Likewise, set-builder notation often uses the same kind of equals-as-iterate.

\displaystyle A = \{ n^2 : n = 1, 2, \dots, 100 \}

In Python, or interpreting the expression literally, the value of n would be a tuple, producing a type error. (In Javascript, it produces 2. How could it be Javascript if it didn’t?)

Next up we have the sloppiness of functions. Let f(x) = 2x + 3. This is a function, and x is a variable. Rather than precisely say, f(2) = 7, we say that for x=2, f(x) = 7. So x is simultaneously an indeterminate input and a concrete value. The same scoping for programming functions bypass the naive expectation that equality means “now and forever.” Couple that with the question-as-equation f(x) = 7, in which one asks what values of x produce this result, if any, and you begin to see how deep the rabbit hole goes. To understand what someone means when they say f(x) = 7, you need to know the context.

But this is just the tip of the iceberg, and we’re drilling deep. The point is that = carries with it all kinds of baggage, not just the scope of a particular binding of a variable.

Continuing with functions, we have rational expressions like f(x) = \frac{(x+1)x}{x}. One often starts by saying “let’s let f be this function.” Then we want to analyze it, and in-so-doing we simplify to f(x) = x+1. To keep ourselves safe, we modify the domain of f to exclude x=0 post-hoc. But the flow of the argument is the same: we defer the exclusion of x=0 until we need it, meaning the equality at the beginning is a different equality than at the end. In effect, we have an infinitude of different kinds of equality for functions, one for each choice of what to exclude from the domain. And a mathematical proof might switch between them as needed.

“Why not just define a new function g with a different domain,” you ask? You can, but mathematicians don’t. And if you’re arguing in favor or against a particular notation, and using “mathematics” as your impenetrable shield, you’ve got to remember the famous definition of Reuben Hersh, that “mathematics is what mathematicians do.” For us, that means you can’t claim superiority based on an idea of mathematics that disagrees with mathematical practice. And mathematics, dear reader, is messier than programmers and philosophers would have one believe.

And now we turn to the Great Equality Contextualizer, the isomorphism. 

You see, all over mathematics there are objects which are not equal, but we want them to be. When you study symmetry, say, you learn that there is an algebraic structure to symmetry called a group. And the same structure—that is, the same true underlying relationships between the symmetries of a thing—can show up in many different guises. As a set, as a picture, as a class of functions, in polynomials and compass constructions and wallpapers, oh my! In each of these things we want to say that two symmetry structures are the same even if they look different. We want to overload equality when four-fold rotational symmetry applies to my table as well as a four-pointed star.

The tool we use for that is called an isomorphism. In brief terms, it’s a function between two objects, with an inverse, that preserves the structure you care about both ways. In fact, there is a special symbol for when two things are isomorphic, and it’s often \cong. But \cong is annoying to write, and it really just means “is the same as” the same way equality does. So mathematicians often drop the squiggle and use =.

Plus, there are a million kinds of isomorphism. Groups, graphs, vector spaces, rings, fields, modules, algebras, rational functions, varieties, Lie groups, *breathe* topological spaces, manifolds of all stripes, sheaves, schemes, lattices, knots, the list just keeps going on and on and on! No way are we making up a symbol for each one of these and the hundreds of variations we might come up with. And moreover, when you say two things are isomorphic, that gives you absolutely no indication of how they are isomorphic. It fact, it can be extremely tedious to compute isomorphisms between things, and it’s even known to be uncomputable in extreme cases! What good is equality if you can’t even check it?

But wait! You might ask, having read this blog for a while and knowing better than to not question a claim. All of these uses of equality are still equivalence relations, and x = x + 1 is not an equivalence relation!

Well, you got me there. Mathematicians love to keep equality as an equivalence relation. When mathematicians need to define an algorithm where the value of x changes in a nontrivial way, it’s usually done by setting x_0 equal to some starting value and letting x_{n} be defined as some function of x_{n-1} and smaller terms, like the good ol’ Fibonacci sequence x_0 = x_1 = 1 and x_n = x_{n-1} + x_{n-2}.

If mutation is so great, why do mathematicians use recursion so much? Huh? Huh?

Well, I’ve got two counterpoints. The first is that the goal here is to reason about the sequence, not to describe it in a way that can be efficiently carried out by a computer. When you say x = x + 1, you’re telling the computer that the old value of x need not linger, and you can do away with the space occupied by the previous value of x. To achieve the same result with recursion requires a whole other can of worms: memoization and tail recursive style and compiler optimizations to shed stack frames. It’s a lot more work to understand all that (to get to an equivalent solution) than it is to understand mutation! Simply stated, the goals of mathematics and programming are quite differently aligned. The former is about understanding a thing, and the latter is more often about describing a concrete process under threat of limited resources.

My second point is that mathematical notation is so flexible and adaptable that it doesn’t need mutation the same way programming languages need it. In mathematics we have no stack overflows, no register limits or page swaps, no limitations on variable names or memory allocation, our brains do the continuation passing for us, and we can rewrite history ad hoc and pile on abstractions as needed to achieve a particular goal. Even when you’re describing an algorithm in mathematics, you get the benefits of mathematical abstractions. A mathematician could easily introduce = as mutation in their work. Nothing stops them from doing so! It’s just that they rarely have a genuine need for it.

But of course, none of this changes that languages could use := or “let” instead of = for assignment. If a strict adherence to asymmetry for asymmetric operations helps you sleep at night, so be it. My point is that the case when = means assignment is an extremely simple bit of context. Much simpler than the albatrossian mental burden required to understand what mathematicians really mean when they write A = B.

Postscript: I hope everyone reading this realizes I’m embellishing a bit for the sake of entertainment. If you want to fight me, tell me the best tree isn’t aspen. I dare you.

Postpostscript: embarrassingly, I completely forgot about Big-O notation and friends (despite mentioning it in the article!) as a case where = does not mean equality! f(n) = O(log n) is a statement about upper bounds, not equality! Thanks to @lreyzin for keeping me honest.

A parlor trick for SET

Tai-Danae Bradley is one of the hosts of PBS Infinite Series, a delightful series of vignettes into fun parts of math. The video below is about the same of SET, a favorite among mathematicians. Specifically, Tai-Danae explains how SET cards lie in (using more technical jargon) a vector space over a finite field, and that valid sets correspond to lines. If you don’t immediately know how this would work, watch the video.

In this post I want to share a parlor trick for SET that I originally heard from Charlotte Chan. It uses the same ideas from the video above, which I’ll only review briefly.

In the game of SET you see a board of cards like the following, and players look for sets.

SetCards

Image source: theboardgamefamily.com

A valid set is a triple of cards where, feature by feature, the characteristics on the cards are either all the same or all different. A valid set above is {one empty blue oval, two solid blue ovals, three shaded blue ovals}. The feature of “fill” is different on all the cards, but the feature of “color” is the same, etc.

In a game of SET, the cards are dealt in order from a shuffled deck, players race to claim sets, removing the set if it’s valid, and three cards are dealt to replace the removed set. Eventually the deck is exhausted and the game is over, and the winner is the player who collected the most sets.

There are a handful of mathematical tricks you can use to help you search for sets faster, but the parlor trick in this post adds a fun variant to the end of the game.

Play the game of SET normally, but when you get down to the last card in the deck, don’t reveal it. Keep searching for sets until everyone agrees no visible sets are left. Then you start the variant: the first player to guess the last un-dealt card in the deck gets a bonus set.

The math comes in when you discover that you don’t need to guess, or remember anything about the game that was just played! A clever stranger could walk into the room at the end of the game and win the bonus point.

Theorem: As long as every player claimed a valid set throughout the game, the information on the remaining board uniquely determines the last (un-dealt) card.

Before we get to the proof, some reminders. Recall that there are four features on a SET card, each of which has three options. Enumerate the options for each feature (e.g., {Squiggle, Oval, Diamond} = {0, 1, 2}).

While we will not need the geometry induced by this, this implies each card is a vector in the vector space \mathbb{F}_3^4, where \mathbb{F}_3 = \mathbb{Z}/3\mathbb{Z} is the finite field of three elements, and the exponent means “dimension 4.” As Tai-Danae points out in the video, each SET is an affine line in this vector space. For example, if this is the enumeration:

joyofset

Source: “The Joy of Set

Then using the enumeration, a set might be given by

\displaystyle \{ (1, 1, 1, 1), (1, 2, 0, 1), (1, 0, 2, 1) \}

The crucial feature for us is that the vector-sum (using the modular field arithmetic on each entry) of the cards in a valid set is the zero vector (0, 0, 0, 0). This is because 1+1+1 = 0, 2+2+2 = 0, and 1+2+3=0 are all true mod 3.

Proof of Theorem. Consider the vector-valued invariant S_t equal to the sum of the remaining cards after t sets have been taken. At the beginning of the game the deck has 81 cards that can be partitioned into valid sets. Because each valid set sums to the zero vector, S_0 = (0, 0, 0, 0). Removing a valid set via normal play does not affect the invariant, because you’re subtracting a set of vectors whose sum is zero. So S_t = 0 for all t.

At the end of the game, the invariant still holds even if there are no valid sets left to claim. Let x be the vector corresponding to the last un-dealt card, and c_1, \dots, c_n be the remaining visible cards. Then x + \sum_{i=1}^n c_i = (0,0,0,0), meaning x = -\sum_{i=1}^n c_i.

\square

I would provide an example, but I want to encourage everyone to play a game of SET and try it out live!

Charlotte, who originally showed me this trick, was quick enough to compute this sum in her head. So were the other math students we played SET with. It’s a bit easier than it seems since you can do the sum feature by feature. Even though I’ve known about this trick for years, I still require a piece of paper and a few minutes.

Because this is Math Intersect Programming, the reader is encouraged to implement this scheme as an exercise, and simulate a game of SET by removing randomly chosen valid sets to verify experimentally that this scheme works.

Until next time!

Earthmover Distance

Problem: Compute distance between points with uncertain locations (given by samples, or differing observations, or clusters).

For example, if I have the following three “points” in the plane, as indicated by their colors, which is closer, blue to green, or blue to red?

example-points.png

It’s not obvious, and there are multiple factors at work: the red points have fewer samples, but we can be more certain about the position; the blue points are less certain, but the closest non-blue point to a blue point is green; and the green points are equally plausibly “close to red” and “close to blue.” The centers of masses of the three sample sets are close to an equilateral triangle. In our example the “points” don’t overlap, but of course they could. And in particular, there should probably be a nonzero distance between two points whose sample sets have the same center of mass, as below. The distance quantifies the uncertainty.

same-centers.png

All this is to say that it’s not obvious how to define a distance measure that is consistent with perceptual ideas of what geometry and distance should be.

Solution (Earthmover distance): Treat each sample set A corresponding to a “point” as a discrete probability distribution, so that each sample x \ in A has probability mass p_x = 1 / |A|. The distance between A and B is the optional solution to the following linear program.

Each x \in A corresponds to a pile of dirt of height p_x, and each y \in B corresponds to a hole of depth p_y. The cost of moving a unit of dirt from x to y is the Euclidean distance d(x, y) between the points (or whatever hipster metric you want to use).

Let z_{x, y} be a real variable corresponding to an amount of dirt to move from x \in A to y \in B, with cost d(x, y). Then the constraints are:

  • Each z_{x, y} \geq 0, so dirt only moves from x to y.
  • Every pile x \in A must vanish, i.e. for each fixed x \in A, \sum_{y \in B} z_{x,y} = p_x.
  • Likewise, every hole y \in B must be completely filled, i.e. \sum_{y \in B} z_{x,y} = p_y.

The objective is to minimize the cost of doing this: \sum_{x, y \in A \times B} d(x, y) z_{x, y}.

In python, using the ortools library (and leaving out a few docstrings and standard import statements, full code on Github):

from ortools.linear_solver import pywraplp

def earthmover_distance(p1, p2):
    dist1 = {x: count / len(p1) for (x, count) in Counter(p1).items()}
    dist2 = {x: count / len(p2) for (x, count) in Counter(p2).items()}
    solver = pywraplp.Solver('earthmover_distance', pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)

    variables = dict()

    # for each pile in dist1, the constraint that says all the dirt must leave this pile
    dirt_leaving_constraints = defaultdict(lambda: 0)

    # for each hole in dist2, the constraint that says this hole must be filled
    dirt_filling_constraints = defaultdict(lambda: 0)

    # the objective
    objective = solver.Objective()
    objective.SetMinimization()

    for (x, dirt_at_x) in dist1.items():
        for (y, capacity_of_y) in dist2.items():
            amount_to_move_x_y = solver.NumVar(0, solver.infinity(), 'z_{%s, %s}' % (x, y))
            variables[(x, y)] = amount_to_move_x_y
            dirt_leaving_constraints[x] += amount_to_move_x_y
            dirt_filling_constraints[y] += amount_to_move_x_y
            objective.SetCoefficient(amount_to_move_x_y, euclidean_distance(x, y))

    for x, linear_combination in dirt_leaving_constraints.items():
        solver.Add(linear_combination == dist1[x])

    for y, linear_combination in dirt_filling_constraints.items():
        solver.Add(linear_combination == dist2[y])

    status = solver.Solve()
    if status not in [solver.OPTIMAL, solver.FEASIBLE]:
        raise Exception('Unable to find feasible solution')

    return objective.Value()

Discussion: I’ve heard about this metric many times as a way to compare probability distributions. For example, it shows up in an influential paper about fairness in machine learning, and a few other CS theory papers related to distribution testing.

One might ask: why not use other measures of dissimilarity for probability distributions (Chi-squared statistic, Kullback-Leibler divergence, etc.)? One answer is that these other measures only give useful information for pairs of distributions with the same support. An example from a talk of Justin Solomon succinctly clarifies what Earthmover distance achieves

Screen Shot 2018-03-03 at 6.11.00 PM.png

Also, why not just model the samples using, say, a normal distribution, and then compute the distance based on the parameters of the distributions? That is possible, and in fact makes for a potentially more efficient technique, but you lose some information by doing this. Ignoring that your data might not be approximately normal (it might have some curvature), with Earthmover distance, you get point-by-point details about how each data point affects the outcome.

This kind of attention to detail can be very important in certain situations. One that I’ve been paying close attention to recently is the problem of studying gerrymandering from a mathematical perspective. Justin Solomon of MIT is a champion of the Earthmover distance (see his fascinating talk here for more, with slides) which is just one topic in a field called “optimal transport.”

This has the potential to be useful in redistricting because of the nature of the redistricting problem. As I wrote previously, discussions of redistricting are chock-full of geometry—or at least geometric-sounding language—and people are very concerned with the apparent “compactness” of a districting plan. But the underlying data used to perform redistricting isn’t very accurate. The people who build the maps don’t have precise data on voting habits, or even locations where people live. Census tracts might not be perfectly aligned, and data can just plain have errors and uncertainty in other respects. So the data that district-map-drawers care about is uncertain much like our point clouds. With a theory of geometry that accounts for uncertainty (and the Earthmover distance is the “distance” part of that), one can come up with more robust, better tools for redistricting.

Solomon’s website has a ton of resources about this, under the names of “optimal transport” and “Wasserstein metric,” and his work extends from computing distances to computing important geometric values like the barycenter, computational advantages like parallelism.

Others in the field have come up with transparency techniques to make it clearer how the Earthmover distance relates to the geometry of the underlying space. This one is particularly fun because the explanations result in a path traveled from the start to the finish, and by setting up the underlying metric in just such a way, you can watch the distribution navigate a maze to get to its target. I like to imagine tiny ants carrying all that dirt.

Screen Shot 2018-03-03 at 6.15.50 PM.png

Finally, work of Shirdhonkar and Jacobs provide approximation algorithms that allow linear-time computation, instead of the worst-case-cubic runtime of a linear solver.

NP-hard does not mean hard

When NP-hardness pops up on the internet, say because some silly blogger wants to write about video games, it’s often tempting to conclude that the problem being proved NP-hard is actually very hard!

“Scientists proved Super Mario is NP-hard? I always knew there was a reason I wasn’t very good at it!” Sorry, these two are unrelated. NP-hardness means hard in a narrow sense this post should hopefully make clear. After that, we’ll explore what “hard” means in a mathematical sense that you can apply beyond NP-hardness to inform your work as a programmer.

When a problem is NP-hard, that simply means that the problem is sufficiently expressive that you can use the problem to express logic. By which I mean boolean formulas using AND, OR, and NOT. In the Super Mario example, the “problem” is a bundle of (1) the controls for the player (2) the allowed tiles and characters that make up a level, and (3) the goal of getting from the start to the end. Logic formulas are encoded in the creation of a level, and solving the problem (completing the level) is the same as finding conditions to make the logical formula true.

mario-clause-gadget

The clause gadget for the original Super Mario Brothers, encoding an OR of three variables.

In this sense, NP-hardness doesn’t make all of Super Mario hard. The levels designed to encode logical formulas are contrived, convoluted, and contorted. They abuse the rules of the game in order to cram boolean logic into it. These are worst case levels. It’s using Mario for a completely unintended purpose, not unlike hacking. And so NP-hardness is a worst case claim.

To reiterate, NP-hardness means that Super Mario has expressive power. So expressive that it can emulate other problems we believe are hard in the worst case. And, because the goal of mathematical “hardness” is to reason about the limitations of algorithms, being able to solve Super Mario in full generality implies you can solve any hard subproblem, no matter how ridiculous the level design.

The P != NP conjecture says that there’s no polynomial time algorithm to determine whether boolean logic formulas are satisfiable, and so as a consequence Super Mario (in full generality) also has no polynomial time algorithm.

That being said, in reality Super Mario levels do not encode logical formulas! If you use the knowledge that real-world Super Mario levels are designed in the way they are (to be solvable, fun), then you can solve Super Mario with algorithms. There are many examples.

In general, the difficulty of a problem for humans is unrelated to the difficulty for algorithms. Consider multiplication of integers. This is a trivial problem for computers to solve, but humans tend to struggle with it. It’s an amazing feat to be able to multiply two 7 digit numbers in less than 5 seconds, whereas computers can multiply two thousand-digit numbers in milliseconds.

Meanwhile, protein folding is known to be an NP-hard problem, but it’s been turned into a game sufficiently easy for humans to solve that players have contributed to scientific research. Indeed, even some of the most typically cited NP-hard problems, like traveling salesman, have heuristic, practical algorithmic solutions that allow one to solve them (very close to optimally) in hours on inputs as large as every city on earth.

So the mathematical notions of hardness are quite disconnected from practical notions of hardness. This is not even to mention that some NP-hard problems can be efficiently approximated to within any desired accuracy.

Let’s dig into the math a bit more. “Hardness” is a family of ideas about comparisons between problems based on reusability of algorithmic solutions. Loosely speaking, a problem R is hard with respect to a class of problems C if an algorithm solving R can be easily transformed into an algorithm solving any problem in C. You have to say what kinds of transformations are allowed, and the transformation can be different for different target problems in C, but that’s the basic idea.

In the Super Mario example, if you want to solve logical formulas, you can transform a hypothetically perfect mario-level-playing algorithm into a logic solver by encoding the formula as a level and running the mario-level-playing algorithm on it as a black box. Add an if statement to the end to translate “level can/can’t be finished” to “formula can/can’t be satisfied,” and the transformation is complete. It’s important for NP-hardness that the transformation only takes polynomial time. Other kinds of hardness might admit more or restrict to fewer resources.

And so this is what makes Mario NP-hard, because boolean logic satisfiability is NP-hard. Any problem in NP can be solved by a boolean logic solver, and hence also by a mario-level-player. The fact that boolean logic solving is NP-hard is a difficult theorem to prove. But if we assume it’s true, you can compose the transformations to get from any NP problem to Super Mario.

As a simple example of a different kind of hardness, you can let C be the class of problems solvable using only a finite amount of memory (independent of the input). You have probably heard of this class of problems by another name, but I’ll keep you guessing until the end of the post. A C-hard problem R is one for which an algorithmic solution can be repurposed to solve any finite-memory-solvable problem.

We have to be careful: if the transformation between solutions allows us polynomial time (in the size of the input) like it did for NP-hardness, then we might have enough time in the transformation alone to solve the entire problem, removing the need for a solution to R in the first place! For this reason, we have to limit the amount of work that can be done in the transformation. We get a choice here that influences how interesting or useful the definition of hardness is, but let’s just pick one and say that the transformation can only use finite time (independent of the input).

To be fair, I actually don’t know if there are any hard problems with respect to this definition. There probably are, but chances are good that they are not members of C, and that’s where the definition of hardness gets really interesting. If you have a problem in C which is also C-hard, it’s called complete for C. And once you’ve found a complete problem, from a theoretical perspective you’re a winner. You’ve found a problem which epitomizes the difficulty of solving problems in C. And so it’s a central aim of researchers studying a complexity class to find complete problems. As they say in the business, “ABC: always be completing.”

As a more concrete and interesting example, the class P of all polynomial-time solvable problems has a complete problem. Here the transformations are a bit up in the air. They could either be logarithmic-space computations, or what’s called NC, which can be thought of as poly-logarithmic time (very fast) parallel computations. I only mention NC because it allows you to say “P-complete problems are hard to parallelize.”

Regardless of the choice, there are a number of very useful problems known to be P-complete. The first is the Circuit Value Problem, given a circuit (described by its gates and wires using any reasonable encoding) and an input to the circuit, what is the output?

Others include linear programming (optimize this linear function with respect to linear constraints), data compression (does the compressed version of a string s using Lempel–Ziv–Welch contain a string t?), and type inference for partial types. There are many more in this compendium of Greenlaw et al. Each one is expressive enough to encode any instance of the other, and any instance of any problem in P. It’s quite curious to think that gzip can solve linear programs, but that’s surely no curiouser than super mario levels encoding boolean logic.

Just as with NP-hardness, when a problem is P-hard that doesn’t automatically mean it’s easy or hard for humans, or that typical instances can’t be easily parallelized. P-hardness is also a worst case guarantee.

Studying P-completeness is helpful in the same way NP-completeness is helpful. Completeness informs you about whether you should hope to find a perfect solution or be content with approximations and heuristics (or incorporate problem context to make it easier). Knowing a problem is P-complete means you should not expect perfect efficient parallel algorithms, or perfect efficient algorithms that use severely limited space. Knowing a problem is NP-hard means you should not expect a perfect polynomial time solution. In other words, if you are forced to work with those restrictions, the game becomes one of tradeoffs. Hardness and completeness focus and expedite your work, and clarify a principled decision making process.

Until next time!

P.S. The class of problems solvable in a finite amount of memory is just the class of regular languages. The “finite memory” is the finite state machine used to solve them.