Those of you who are familiar with Garey and Johnson’s 1979 classic, Computers and Intractability: a guide to the theory of NP-completeness, may notice I’m simply “porting” their introduction, including the dialogue, to the statistics world.
Imagine Andrew had tasked me and Matt Hoffman with fitting simple standard (aka isotropic, aka spherical) Gaussian mixtures rather than hierarchical models. Let’s say that Andrew didn’t like that K-means got a different answer every time he ran it, K-means++ wasn’t much better, and even using soft-clustering (i.e., fitting the stat model with EM) didn’t let him replicate simulated data. Would we have something like Stan for mixtures. Sadly, no. Matt and I may have tried and failed. We wouldn’t want to go back to Andrew and say,
- “We can’t find an efficient algorithm. I guess we’re just too dumb.”
We’re computer scientists and we know about proving hardness. We’d like to say,
- “We can’t find an efficient algorithm, because no such algorithm is possible.”
But that would’ve been beyond Matt’s and my grasp, because, in this particular case, it would require solving the biggest open problem in theoretical computer science. Instead, it’s almost certain we would have come back and said,
- “We can’t find an efficient algorithm, but neither can all these famous people.”
That seems weak. Why would we say that? Because we could’ve proven that the problem is NP-hard. A problem is in the class P if it can be solved in polynomial time with a deterministic algorithm. A problem is in the class NP when there is a non-deterministic (i.e., infinitely parallel) algorithm to solve it in polynomial time. It’s NP-hard if it’s just as hard as any other NP algorithm (formally specified through reductions, a powerful CS proof technique that’s the basis of Gödel’s incompleteness theorem). An NP-hard algorithm often has a non-deterministic algorithm to solve it makes a complete set of (exponentially many) guesses in parallel and then spends polynomial time on each one verifying whether or not it is a solution. An algorithm is NP-complete if it is NP-hard and a member of NP. Some well known NP-complete problems are bin packing, satisfiability in propositional logic, and the traveling salesman problem—there’s a big list of NP-complete problems.
Nobody has found a tractable algorithm to solve an NP-hard problem. When we (computer scientists) say “tractable,” we mean solvable in polynomial time with a deterministic algorithm (i.e., the problem is in P). The only known algorithms for NP-hard problems are exponential. Researchers have been working for the last 50+ years trying to prove that the class of NP problems is disjoint from the class of P problems.
In other words, there’s a Turing Award waiting for you if you can actually turn response (3) into response (2).
In the case of mixtures of standard (spherical, isotropic) Gaussians there’s a short JMLR paper with a proof that maximum likelihood estimation is NP-hard.
- Tosh and Dasgupta. 2016. Maximum Likelihood Estimation for Mixtures of Spherical Gaussians is NP-hard. JMLR
And yes, that’s the same Tosh as who was the first author of the “piranha” paper.
Ising models that are not restricted to be planar are also NP-hard.
- Lucas. 2014. Ising formulations of many NP problems. Frontiers in Physics.
What both these problems have in common is that they are combinatorial and require inference over sets. I think (though am really not sure) that one of the appeals of quantum computing is potentially solving NP-hard problems.
P.S. How this story really would’ve went is that we would’ve told Andrew that some simple distributions over NP-hard problem instances lead to expected polynomial time algorithms and we’d be knee-deep in the kinds of heuristics used to pack container ships efficiently.
Worst-case complexity is NP-hard, but for K-means I’m not sure that’s actually a useful metric. For example, if you add a little bit of gaussian noise to the data you’re down in polynomial time: https://arxiv.org/abs/0904.1113.
As the saying goes, uncertainty greases the wheels of commerce!
As far as I can tell (the language they use is pretty unfamiliar to me), it’s just saying that the run time of the algorithm to find a local mode is polynomial. I couldn’t find anything that claims to be able to find the MLE (the global maximum).
I was specifically talking about the problem where we assume each cluster has an isotropic Gaussian distribution and we want to find the MLE. K-means is an EM algorithm for this problem that is attenuated. You can run a regular EM algorithm, but that also won’t find the MLE for you in hard problem instances.
I had thought K-means was NP-complete (i.e., NP-hard and in the class NP). See, for example, this paper: https://cseweb.ucsd.edu/~avattani/papers/kmeans_hardness.pdf
In that case, proving that K-means is exponentially hard would settle P != NP, so that can’t be what’s going on.
I resonate with this.
Around 1970, I was starting on PhD and advisor suggested developing an optimal algorithm for a 2-machine, n jobs,open shop scheduling problem.
I wrote a lot of APL to test out different algorithms, BUT:
Of course, this was just before NP started to be understood, and it turns out that the general problem is NP-hard, which people proved a few years later.
https://en.wikipedia.org/wiki/Open-shop_scheduling
Luckily, there were other things to do.
I vaguely recall hearing of other students at that time who started such projects only to find they were NP.
NP-hardness turns out to be an unsatisfying measure of complexity because it’s only looking at the hardest instances of a problem, and those often turn out to be a lot harder than typical instances of the problem. There’s been a lot of work in the past couple decades on coming up with more refined notions of complexity, and there’s a recent book (https://www.cambridge.org/core/books/beyond-the-worstcase-analysis-of-algorithms/8A8128BBF7FC2857471E9CA52E69AC21) that surveys the field. The last section of the book on applications to ML and statistics may be of particular interest to readers of this blog.
Right. That’s what I was getting at when I wrote “some simple distributions over NP-hard problem instances lead to expected polynomial time algorithms” in the last sentence of the post. You have to make an assumption about the probability of problem instances to make claims about expected running time. And you still don’t wind up with a general algorithm that can solve any K-means problem in polynomial time.
If we’re talking about worst case complexity (as other commenters have pointed out, this isn’t always the right thing to think about in practice), then my understanding is that quantum computers are not likely to help with NP-hard problems (in the sense of having quantum algorithms that run in polynomial time). Grover’s algorithm can give a quadratic speedup for NP-hard problems like SAT (which could eventually matter if we ever make quantum computing practical), but does not lead to a quantum polynomial time algorithm.
The famous result of giving exponential speedup over classical computing (in theory) is Shor’s algorithm which can be used to solve the factoring problem in quantum polynomial time. The issue here is that factoring is not known to be NP-hard, and complexity theorists believe this is likely not the case.
Larry Wasserman’s infamous remark, “I have decided that mixtures, like tequila, are inherently evil and should be avoided at all costs,” is entirely justified as a general principle. But some infinite mixtures, e.g. Gaussian location mixtures, have unique, easily computed NPMLEs as shown here: https://www.tandfonline.com/doi/abs/10.1080/01621459.2013.869224, and in earlier work of Bruce Lindsay.
Thanks for that quote. I feel the same way about tequila :-). I’m so used to working in general black-box settings that I often forget you can actually make a lot of progress on special cases.
I’m surprised no-one mentioned how MLE is basically trash and I’d take 10 iid samples from a posterior distribution over an MLE in a heartbeat, and that getting MCMC samples from a mixture posterior is really not that hard.
Unfortunately, nobody knows how to take 10 i.i.d. draws from a high-dimensional Gaussian mixture conditioned on data! It’s just as intractable as the MLE. Every sampling method I know will get stuck in local modes or at most explore a handful of a prohibitively combinatorial number of modes. Some people, including Yuling Yao in his stacking paper with Andrew, have argued that if you run a whole bunch of MCMC from different starting points and stack the results, it’s better than just running MCMC. I can believe that, but it’s still not truly sampling from the posterior.
Bob, it’s true it’s not easy (ie. impossible in general) to get 10 iid samples, but for tons of models it’s not that hard to get the informational equivalent of 10 IID samples out of well specified mixture models.
Pulling a sample with 10 ESS from say 1/3 * normal(a,b) + 2/3 * normal(c,d) with let’s say a constrained to be near 0, b constrained to be near 10, c constrained to be near 1 and d constrained to be near 2 and let’s say 10 data points drawn from the actual model is pretty easy for example.
Where things usually go wrong is when people have label-switching type symmetries or similar issues. But those are usually model-specification related issues. Sure, sometimes we don’t know how to specify a better model, but now we’re talking about model specification rather than model fitting.
In the context of a “classical” or “likelihoodist” fit, they’re relying on the maximum likelihood point with *no prior* to fit the model. In the above specified model with no priors you just say “there are two gaussians and they have weights q and (1-q) for q between 0 and 1… you can’t fit it, it’s not enough information. But with the priors the whole thing works fine.
In an unconstrained MLE type fit what you often find is that the landscape is very flat in some big region. This is a huge problem for “maximum likelihood” type fitting. The parameter vector can move a long way and produce almost no change in the likelhood, the gradient is flat, finding the single maximum point is hard. it’s not a problem at all for sampling, you get samples inside that big region, and find out that a lot of different models could work.
Yes, when you have many isolated bumps in multi-dimensional space you can have a big problem. But this is usually a model specification issue rather than a sampling issue. Like with “label switching” type symmetries.
Basically what I’m saying is that if you expect “the machinery” to give you “the one true estimate” you have lost the thread before you’ve even started. It’s the wrong problem to solve.