The Remarkable k-means++

1. The Problem

One of the most popular algorithms for clustering is k-means. We start with data vectors . We choose vectors — cluster centers — to minimize the error

Unfortunately, finding to minimize is NP-hard. The usual iterative method, \hrefnosnap{http://en.wikipedia.org/wiki/Lloyd is easy to implement but it is unlikely to come close to minimizing the objective function. So finding

isn’t feasible.

To deal with this, many people choose random starting values, run the -means clustering algorithm then rinse, lather and repeat. In general, this may work poorly and there is no theoretical guarantee of getting close to the minimum. Finding a practical method for approximately minimizing is thus an important practical problem.

2. The Solution

David Arthur and Sergei Vassilvitskii came up with a wonderful solution in 2007 known as k-means++.

The algorithm is simple and comes with a precise theoretical guarantee.

The first step is to choose a data point at random. Call this point . Next, compute the squared distances

Now choose a second point from the data. The probability of choosing is . Now recompute the distance as

Now choose a third point from the data where the probability of choosing is . We continue until we have points . Finally, we run -means clustering using as starting values. Call the resulting centers .

Arthur and Vassilvitskii prove that

The expected value is over the randomness in the algorithm.

There are various improvements to the algorithm, both in terms of computation and in terms of getting a sharper performance bound.

This is quite remarkable. One simple fix, and an intractable problem has become tractable. And the method comes armed with a theorem.

3. Questions

1. Is there an R implementation? It is easy enough to code the algorithm but it really should be part of the basic k-means function in R.
2. Is there a version for mixture models? If not, it seems like a paper waiting to be written.
3. Are there other intractable statistical problems that can be solved using simple randomized algorithms with provable guarantees? (MCMC doesn’t count because there is no finite sample guarantee.)

4. Reference

Arthur, D. and Vassilvitskii, S. (2007). k-means++: The advantages of careful seeding. Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms. 1027–1035.

Advertisement

Machine Learning In New York

I wanted to let you know that, on October 19, the New York Academy of Sciences is hosting its
7th Annual Machine Learning Symposium.

For details: see here:

link to symposium

Screening: Everything Old is New Again

Screening: Everything Old Is New Again

Screening is one of the oldest methods for variable selection. It refers to doing a bunch of marginal (single covariate) regressions instead of one multiple regression. When I was in school, we were taught that it was a bad thing to do.

Now, screening is back in fashion. It’s a whole industry. And before I throw stones, let me admit my own guilt: see Wasserman and Roeder (2009).

1. What Is it?

Suppose that the data are with

To simplify matters, assume that , and . Let us assume that we are in the high dimensional case where . To perform variable selection, we might use something like the lasso.

But if we use screening, we instead do the following. We regress on , then we regress on , then we regress on . In other words, we do one-dimensional regressions. Denote the regression coefficients by . We keep the covariates associated with the largest values of . We then might do a second step such as running the lasso on the covariates that we kept.

What are we actually estimating when we regression on the covariate? It is easy to see that

where

and is the correlation between and .

2. Arguments in Favor of Screening

If you miss an important variable during the screening phase you are in trouble. This will happen if is big but is small. Can this happen?

Sure. You can certainly find values of the ‘s and the to make big and make small. In fact, you can make huge while making . This is sometimes called unfaithfulness in the literature on graphical models.

However, set of vectors that are unfaithful has Lebesgue measure 0. Thus, in some sense, unfaithfulness is “unlikely” and so screening is safe.

3. Arguments Against Screening

Not so fast. In order to screw up, it is not necessary to have exact unfaithfulness. All we need is approximate unfaithfulness. And the set of approximately unfaithful ‘s is a non-trivial subset of .

But it’s worse than that. Cautious statisticians want procedures that have properties that hold uniformly over the parameter space. Screening cannot be successful in any uniform sense because of the unfaithful (and nearly unfaithful) distributions.

And if we admit that the linear model is surely wrong, then things get even worse.

4. Conclusion

Screening is appealing because it is fast, easy and scalable. But it makes a strong (and unverifiable) assumption that you are not unlucky and have not encountered a case where is small but is big.

Sometimes I find the arguments in favor of screening to be appealing but when I’m in a more skeptical (sane?) frame of mind, I find screening to be quite unreasonable.

What do you think?

Wasserman, L. and Roeder, K. (2009). High dimensional variable selection. Annals of statistics, 37, 2178.

High Dimensional Undirected Graphical Models

High Dimensional Undirected Graphical Models
Larry Wasserman

Graphical models have now become a common tool in applied statistics. Here is a graph representing stock data:

Here is one for proteins: (Maslov and Sneppen, 2002).

And here is one for the voting records of senators: (Banerjee, El Ghaoui, and d’Aspremont ,2008).

Like all statistical methods, estimating graphs is challenging in high dimensional problems.

1. Undirected Graphical Models

Let be a random vector with distribution . A graph for has nodes, or vertices, one for each variable. Some pairs of vertices are connected by edges. We can represent the edges as a set of unordered pairs: if and only if there is an edge between and .

We omit an edge between and to mean that and are independent, given the other variables which we write as

where “rest” denotes the rest of the variables. For example, in this graph

we see that since there is not edge between and .

Now suppose we observe random vectors

The goal is to estimate from the data. As an example, imagine that is a vector of gene expression levels for subject . The graph gives us an idea of how the genes are related.

2. The Gaussian Case

Things are easiest if we assume that is a multivariate Gaussian. Suppose . In this case, there is no edge between and iff where . If the dimension of is smaller than the sample size , then estimating the graph is straightforward. We can use the sample covariance matrix

to estimate and we use to estimate . It is then easy to test . We put an edge between and if the test rejects. The R package SIN will do this for you.

When , the simple approach above won’t work. Perhaps the easiest approach is the method due to Meinshausen and B{ü}hlmann (2006) which John Lafferty has aptly named parallel regression. The idea is this: you regress on all the other variables, then you regress on all the other variables, etc. For each regression, you use the lasso which yields sparse estimators. When you regress on the others, there will be a regression coefficient for . Call this . Similarly, when you regress on the others, there will be a regression coefficient for . Call this . If and then you put an edge between and .

An alternative — the glasso (graphical lasso) — is to maximize the Gaussian log-likelihood with a penalty:

The resulting estimator is sparse: many of its elements are 0. The non-zeroes denote edges. A fast implementation in R due to Han Liu can be found here.

3. Relaxing Gaussianity

The biggest drawback of the glasso is the assumption of Gaussianity. With my colleagues John Lafferty and Han Liu and others, I have done some work on more nonparametric approaches.

Let me briefly describe the approach in Liu, Xu, Gu, Gupta, Lafferty and Wasserman (2011). The idea is to restrict the graph to be a forest, which is a graph with no cycles.

When is a forest, the density can be written as

where is the set of edges. We can estimate the density by inserting estimates of the univariate and bivariate marginals:

But to find the graph we need to estimate the edge set .

Any forest defines a density by the above equation. If the true density were known, we could choose to minimize the Kullback-Leibler distance

This maximization can be done by Kruskal’s algorithm (Kruskal 1956) also known, in this context, as the Chow-Liu algorithm (Chow and Liu 1968). It works like this.

For each pair let

be the mutual information between and . We start with an empty tree and add edges greedily according to the value of . First we connect the pair of variables with the largest and then the second largest and so on. However, we do not add an edge if it forms a cycle.

Of course, we don’t know so, instead, we use the data to estimate the mutual information . If we simply run the Chow-Liu algorithm, we get a tree, that is, a fully connected forest. But we also use a hold-out sample to decide on when to stop adding edges. And that’s it.

There are other ways to relax the Gaussian assumption. For example, here is an approach that uses rank statistics and copulas.

4. The Future?

Not too long ago, estimating a high dimensional graph was unthinkable. Now they are used routinely. The biggest thing missing is a good measure of uncertainty. One can do some obvious things, such as bootstrapping, but it is not clear that the output is meaningful.

5. References

Banerjee, O. and El Ghaoui, L. and d’Aspremont, A. (2008). Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data. The Journal of Machine Learning Research, 9, 485-516.

Chow, C. and Liu, C. (1968). IEEE Transactions on Information Theory, 14, 462-467.

Kruskal, J.B. (1956). On the shortest spanning subtree of a graph and the traveling salesman problem. Proceedings of the American Mathematical society, 7, 48-50.

Liu, H., Xu, M., Gu, H., Gupta, A., Lafferty, J. and Wasserman, L. (2011). Forest Density Estimation. Journal of Machine Learning Research, 12, 907-951.

Maslov, S. and Sneppen, K. (2002). Specificity and stability in topology of protein networks. Science, 296, 910.

Meinshausen, N. and B{ü}hlmann, P. (2006). High-dimensional graphs and variable selection with the lasso. The Annals of Statistics, 34, 1436-1462.

Hunting for Manifolds

Hunting For Manifolds
Larry Wasserman

In this post I’ll describe one aspect of the problem of estimating a manifold, or manifold learning, as it is called.

Suppose we have data and suppose that the data lie on, or close to, a manifold of dimension . There are two different goals we might have in mind.

The first goal is to use the fact that the data lie near the manifold as a way of doing dimension reduction. The popular methods for solving this problem include isomap, local linear embedding, diffusion maps, and Laplacian eigenmaps among others.

The second goal is to actually estimate the manifold .

I have had the good fortune to work with two groups of people on this problem. With the first group (Chris Genovese, Marco Perone-Pacifico, Isa Verdinelli and me) we have focused mostly on certain geometric aspects of the problem. With the second group (Sivaraman Balakrishnan, Don Sheehy, Aarti Singh, Alessandro Rinaldo and me) we have focused on topology.

It is this geometric problem (my work with the first group) that I will discuss here. In particular, I will focus on these two papers: here and here.

1. Statistical Model

To make some progress, we introduce a statistical model for the data. Here are two models.

Model I (Clutter Noise): With probability we draw from a distribution supported on . With probability we draw from a uniform distribution on some compact set . To be concrete, let’s take . The distribution of is

where is uniform on .

Model II (Additive Noise): We draw from a distribution supported on . Then we let where is a -dimensional Gaussian distribution. In this case, has density given by

where is a Gaussian. We can also write this as

where denotes convolution.

There are some interesting complications. First, note that is a singular distribution. That is, there are sets with positive probability but 0 Lebesgue measure. For example, suppose that is a circle in and that is uniform on the circle. Then the circle has Lebesgue measure 0 but . Model II has the complication that it is a convolution and convolutions are not easy to deal with.

Our goal is to construct an estimate .

2. Hausdorff Distance

We will compare to using the Hausdorff distance. If and are sets, the Hausdorff distance between and is

where

and is a ball of radius centered at .

Let’s decode what this means. Imagine we start to grow the set by placing a ball of size around each point in . We keep growing until it engulfs . Let be the size of needed to do this. Now grow until it engulfs . Let be the size of needed to do this. Then . Here is an example:

There is another way to think of the Hausdorff distance. If is a set and is a point, let be the distance from to :

We call the distance function. Then

Using as our loss function, our goal is to find the minimax risk:

where the infimum is over all estimators and the supremum is over a set of distributions . Actually, we will consider a large set of manifolds and we will associate a distribution with each manifold . This will define the set .

3. Reach

We will consider all -dimensional manifolds with positive reach. The reach is a property of a manifold that appears in many places in manifold estimation theory. I believe it first appeared in Federer (1959).

The reach is is the largest number such that each point in has a unique projection on . A set of non-zero reach has two nice properties. First, it is smooth and second, it is self-avoiding.

Consider a circle in the plane with radius . Every point on the plane has a unique projection on with one exception: the center of the circle is distance from each point on . There is no unique closest point on to . So the reach is .

A plane has infinite reach. A line with corner has 0 reach. Here is a example of a curve in two-dimensions. Suppose you put a line, perpendicular to , at each point on . If the length of the lines is less than reach(M), then the lines won’t cross. If the length of the lines is larger than reach(M), then the lines will cross.

Requiring to have positive reach is a strong condition. There are ways to weaken the condition based on generalization so reach, but we won’t go into those here.

The set of manifolds we are interested in is the set of -manifolds contained in with reach at least . Here, is some arbitrary, fixed positive number.

4. The Results

Recall that that the minimax risk is

In the case of clutter noise it turns out that

Note that the rate depends on but not on the ambient dimension . But in the case of additive Gaussian noise,

The latter result says that approaches 0 very slowly. For all practical purposes, it is not possible to estimate . This might be surprising. But, in fact, the result is not unexpected.

We can explain this by analogy with another problem. Suppose you want to estimate from observations . if satisfies standard smoothness assumptions, then the minimax risk has the form

which is a nice, healthy, polynomial rate. But if we add some Gaussian noise to each observation, the minimax rate of convergence for estimating is logarithmic (Fan 1991). The same thing happens if we add Gaussian noise to the covariates in a regression problem. A little bit of Gaussian noise makes these problems essentially impossible.

Is there a way out of this? Yes. Instead of estimating we estimate a set that is close to . In other words, we have to live with some bias. We can estimate , which we call a surrogate for , with a polynomial rate.

We are finishing a paper on this idea and I’ll write a post about it when we are done.

5. Conclusion

Estimating a manifold is a difficult statistical problem. With additive Gaussian noise, it is nearly impossible. And notice that I assumed that the dimension of the manifold was known. Things are even harder when is not known.

However, none of these problems are insurmountable as long as we chnage the goal from “estimating ” to “estimating an approximation of ” which, in practice, is usually quite reasonable.

6. References

Fan, J. (1991). On the optimal rates of convergence for nonparametric deconvolution problems. Ann. Statist. 19 1257-1272.

Federer, H. (1959). Curvature measures. Trans. Amer. Math. Soc. 93 418-491.

Genovese, C., Perone-Pacifico, M., Verdinelli, I. and Wasserman, L. (2012). Minimax Manifold Estimation. Journal of Machine Learning Research, 13, 1263–1291.

Genovese, C., Perone-Pacifico, M., Verdinelli, I. and Wasserman, L. (2012). Manifold estimation and singular deconvolution under Hausdorff loss. The Annals of Statistics, 49, 941-963.

Robins and Wasserman Respond to a Nobel Prize Winner Continued: A Counterexample to Bayesian Inference?

Robins and Wasserman Respond to a Nobel Prize Winner Continued: A Counterexample to Bayesian Inference?

This is a response to Chris Sims’ comments on our previous blog post. Because of the length of our response, we are making this a new post rather than putting it in the comments of the last post.

Recall that we observe iid observations , where and are Bernoulli and independent given .
Define and . We assume that is a known function. Also the marginal
density of (with ) is known and uniform on the unit cube in . Our goal is estimation of

The likelihood

factors into two parts – the first depending on and the second on .

1. Selection Bias

This is a point of agreement. There is selection bias if and only if

Note that

Hence, if then
the sample average of in the subset whose is observed (R=1) is unbiased for

In this case, inference is easy for the Bayesian and the frequentist and there is no issue. So we all agree that the interesting case is where there is selection bias,
that is, where

2. Posterior Dependence on

If the prior on the functions and is such that then the posterior does not depend on and the posterior for will not concentrate around the true value of . Again, we believe we all agree on this point.

We note that no one, Bayesian or frequentist, has ever proposed using an estimator that does not depend on in the selection bias case, i.e when is non-zero. (See addendum for more on this point.)

3. Prior Independence Versus Selection Bias

Reading Chris’ comments, the reader might get the impression that prior independence rules out selection bias, i.e. that is,

Therefore, one might conclude that if we want to discuss the interesting case where there is selection bias, then we cannot have .

But this is incorrect. does not imply that To see this, consider the following example.

Suppose that is one dimensional and a Bayesian’s prior for depends only on the two parameters as follows:

where is uniform on and is uniform on .
Then, clearly and are independent under . However, recalling is uniform so we have that for for any fixed ,

Hence

since and are both positive with probability 1.

4. Other Justifications For Prior Dependence?

Since prior independence of and does not imply “no selection bias,” one might instead argue that it is practically unrealistic to have . But we now show that it is realistic.

Suppose a new HMO needs to estimate the fraction of its patient population that will have a MI in the next year, so as to determine the number of cardiac unit beds needed. Each HMO member has had 300 potential risk factors measured: age, weight height, blood pressure, multiple tests of liver, renal, pulmonary, and cardiac function, good and bad cholesterol, packs per day smoked, years smoked, etc. (We will get to 100,000 once routine genomic testing becomes feasible). A general epidemiologist had earlier studied risk factors for MI
by following 5000 of the 50,000 HMO members for a year. Because MI is a rare event, he oversampled subjects whose , in his opinion, indicated a
smaller probability of an MI (. Hence the
sampling fraction was a known, but complex function chosen so as to try to make and negatively correlated.

The world’s leading heart expert, our Bayesian, was hired to estimate based on distribution of in HMO members and the data from the study.
As world’s expert, his beliefs about the risk function would not change upon learning as only reflects a nonexpert’s beliefs. Hence and are a priori independent. Nonetheless, knowing that the epidemiologist had carefully read the expert literature on risk factors for MI, he also believes with high probability that epidemiologist succeeded in having the random variables and be negatively correlated.

What’s more, Robins and Ritov (1997) showed that, if before seeing the data, any Bayesian, cardiac expert or not, thoroughly queries the epidemiologist
(who selected ) about the epidemiologist’s reasoned opinions concerning (but not about ), the Bayesian will then have independent priors. The idea is that once you are satisfied that you have learned from the epidemiologist all he knows about that you did not, you will have an updated prior for
. Your prior for now updated) cannot then change if you subsequently are told Hence, we could take as many Bayesians as you please and arrange it so all had and apriori independent. This last argument is quite general, applying to many settings.

5. Alternative Interpretation

An alternative reading of Chris’s third response and his subsequent post is that, rather than placing a joint prior over the functions
and as above,
his prior is placed over the joint distribution of the random variables and .
If so, he is then correct that making and independent with prior probability one
also implies and thus no selection bias.

However, it appears that from this, he concludes that selection bias, in itself, licenses the dependence of his posterior on .
This is incorrect. As noted above, it is prior dependence of and that licenses posterior dependence on – not prior dependence of and
. Were he correct, our Bayesian cardiac expert’s prior on could have changed upon learning the epidemiologists .

6. What If We Do Use a Prior That Depends on ?

In the above scenario, should ot depend on . But suppose, for whatever reason, one insists on letting depend on .

That still does not mean the posterior will concentrate. Having an estimator that depends on is necessary, but not sufficient, to get consistency and fast rates. It is not enough to use a prior that is a function of . The prior still has to be carefully engineered to ensure that the posterior for will concentrate around the truth.

Chris hints that he can construct such a prior but does not provide an explicit algorithm nor an argument as to why the estimator would be expected to be locally semiparametric efficient. However, it is simple to construct a -consistent
locally semiparametric efficient Bayes estimator as follows.

We tentatively model as a finite dimensional parametric function
with either a smooth or noninformative prior on the parameters , where we take

, and the are basis functions. Then the posterior mean
of will have the same asymptotic distribution as the locally semiparametric efficient regression estimator of Scharfstein et al. (1999) described in our original post. Note that the estimator is consistent, even if the model is wrong.

Of course, this estimator is a clear case of frequentist pursuit Bayes.

7. Conclusion

Here are the main points:

1. If then the posterior will not concentrate.
   Thus, if a Bayesian wants the posterior for to concentrate around the true value,
   he must justify having a prior that is a function of .

2. does not imply an absence of selection bias.
   Therefore, an argument of the form: “we want selection bias so we cannot have prior independence” fails.

3. One can try to argue that prior independence is unrealistic. But as we have shown, this is not the case.
4. But, if after all this, we do insist on letting depend on ,
   it is still not enough. Dependence on is necessary but not sufficient.

We conclude Bayes fails in our example unless one uses a special prior designed just to mimic the frequentist estimator.

8. Addendum: What happens If The Estimator Does Note Depend on ?

The theorem of Robins and Ritov, quoted in our initial post, says that no uniformly consistent estimator that does not depend on can exist in the model which contains all measurable and subject
to with probability 1. Take for
concreteness. In fact, even when we assume
and are quite smooth, there will be little
improvement in performance.

Given has 100,000 dimensions, we can ask how many derivatives and must
and have so that it is possible to construct an estimator of , not depending on that
converges at rate uniformly to over a submodel . Robins et al. (2008) show that it is necessary and sufficient that + and provide an
explicit estimator. More generally, if + derivatives with , the optimal rate is which is approximately when is small compared to An explicit estimator
is constructed in Robins et al (2008) ; Robins et al (2009) prove that the rate cannot be improved on. Given these asymptotic mathematical results, we
doubt any reader can exhibit an estimator, not depending on that will have reasonable finite sample performance under
model or even with, say,
and a sample size of 5,000. By reasonable finite sample performance, we mean an interval estimator that will cover the true at least 95% of the time and that has average length less than or equal to intervals estimators
centered on the improved HT estimators. Nonetheless, we
await any candidate estimators, accompanied by at least some simulation
evidence backing up your claim.

9. References

1. Robins JM, Tchetgen E, Li L, van der Vaart A. (2009).
   Semiparametric Minimax Rates. Electron. J. Statist. Volume 3 (2009),
   1305-1321.

2. Robins JM, Li L, Tchetgen E, van der Vaart A. (2008). Higher order influence
   functions and minimax estimation of nonlinear functionals. Probability and
   Statistics: Essays in Honor of David A. Freedman 2:335-421

3. Robins JM, Ritov Y. (1997). Toward a curse of dimensionality appropriate
   (CODA) asymptotic theory for semi-parametric models. Statistics in Medicine,
   16:285-319.