Why can't one generalize the Kolmogorov-Smirnov test to 2 or more dimensions? The question says it all. I've read both that one can't generalize KS to a dimension equal or larger than two, and that famous implementations like that in Numerical Recipes are simply wrong.
Could you please explain why is so?
 A: I believe it's legitimate to quote the relevant portion of the paragraph in question:

3.  The KS test can not be applied in two or more dimensions. Astronomers often have datasets with points distributed in a plane or higher dimensions, rather than along a line.  Several papers in the astronomical literature purport to present a two-dimensional KS test, and one is reproduced in the famous volume Numerical Recipes.  However, no EDF-based test (this includes KS, AD and related tests) can be applied in two or higher dimensions, because there is no unique way to order the points so that distances between well-defined EDFs can be computed.  One can construct a statistic based on some ordering procedure, and then compute the supremum distances between two datasets (or one dataset and a curve).  But the critical values of the resulting statistic are not distribution-free.  

As stated, this seems too strong.
1) The bivariate distribution function, which is $F(x_1,x_2) = P(X_1\leq x_1,X_2\leq x_2)$ is a map from $\mathbb{R}^2$ to $[0,1]$. That is, the function takes univariate real values between 0 and 1. Those values - being probabilities - are certainly "ordered" already - and this (the value of the function) is the thing we need to make comparisons on for ECDF-based tests. Similarly, the ecdf, $\hat F$ is perfectly well defined in the bivariate case.
I don't think there's necessarily a need to try to turn it into some function of a univariate combined variable as the text suggests. You simply compute $F$ and $\hat F$ at every required combination and compute the difference. 
2) However, on the question of whether it's distribution-free, they have a point: 
a) clearly such a test statistic would not be altered by changes to transformations of the margins, which is to say, if constructed as a test of bivariate independent uniforms, $\mathbf{U}=(U_1,U_2)$, then it works equally well as a test of independent $(X_1,X_2)$ where $U_i=F_i(X_i)$. In that sense, it's distribution-free (we might say 'margin-free').
b) however, there's an underlying point more generally in the broader sense that a naive version of the KS statistic (such as I just described) is not more generally distribution free; we can't simply transform $U$ arbitrarily $X^* = \mathbf{g}(\mathbf{U})$.
In an earlier version of my answer I said: 

There's no difficulty, no problem 

That's wrong. There are indeed issues if there's a change not just of the margins from bivariate independent uniforms, as just mentioned. However, those difficulties have been considered in several ways in a number of papers that yield bivariate/multivariate versions of Kolmogorov-Smirnov statistics that don't suffer from that problem.
I may come back and add some of those references and some discussion of how they work as soon as time permits.
