Concentration of sum of pairwise squared Euclidean distances of random vectors Let $X_1, \ldots, X_n $ be independent random vectors in $B(0, D) \subset R^d$ ($\ell_2$ ball of radius $D$ centered at the origin). I am trying to find the concentration of the following quantity around its expectations:
$$f(X_1, \ldots, X_n) = \frac{1}{n^2} \sum_{1 \leq i,j \leq n} \|X_i - X_j \|_2^2$$
Using McDiarmid's inequality, I can show that 
$$ \Pr(|f(X_1,\ldots,X_n) - \mathbb{E} f(X_1,\ldots,X_n)| \geq \epsilon) \leq 2 \exp\left(- \frac{2n\epsilon^2}{16D^4}\right). $$
This tells me that with high probability, $\epsilon \sim O(1/\sqrt{n})$. 
Usually sum of $n$ independent random variables concentrate with the similar dependence on $n$ (that is, $O(1/\sqrt{n})$) unless I use Berstein's inequality. And it is believed that in some sense, among functions of independent random variables, sums are the least concentrated.
In my situation, the function is a sum of $n^2$ non-independent random variables (if you think of $Z_{ij} = \| X_i - X_j \|_2^2$ as a random variable). So I believe that the concentration should be better than $O(1/\sqrt{n})$.
Does anybody have any idea of how such a guarantee can be achieved? Or can anybody show that the $O(1/\sqrt{n})$ is the best I can hope for? 
 A: You're right, you can get a better concentration inequality with $\epsilon \sim O(1/n)$ considering that's a sum of $n^2$ dependent random variables. You first need to use decoupling for $U$-statistics that relates the dependent setting to the independent one as it is shown in  [1]. In your case, you can show that there exist a constant $C$ such that
$ P\left(\left| \frac{1}{n^2} \sum_{1 \leq i \ne j \leq n} \|X_i - X_j \|_2^2
 - \mathbb{E}f\right| \ge t\right) \le \\
\qquad \qquad C P\left(C\left|\frac{1}{n^2} \sum_{1 \leq i \ne j \leq n} \|X_i^{(1)} - X_j^{(2)} \|_2^2
 - \mathbb{E}f\right| \ge t\right) , $
where $X_i^{(1)}, X_i^{(2)}$ are independent copies of $X_i$.
Then you can apply the McDiarmid inequality on the sum of $n(n-1)$ independent terms and get an inequality achieving the rate $O(1/n)$.
Moreover, as your U-statistic is symmetric you can reverse this inequality to get a lower bound and convince you that you can't do better.
[1] de la Peña, V. H., & Montgomery-Smith, S. J. (1995). Decoupling inequalities for the tail probabilities of multivariate U-statistics. The Annals of Probability, 806-816. Avialable : http://arxiv.org/pdf/math/9309211v2.pdf
A: I wanted to comment but don't have enough reputation yet. With regards to showing the optimality of the $O(1/\sqrt{n})$ rate, can you try the following strategy? Fix $d = 1$ and $D = 1$. Let $X_1, \ldots, X_n$ be iid Rademacher random variables (i.e taking values $1$ and $-1$ each with probability $1/2$). Then the random variable $f(X_1, \ldots, X_n)$ simplifies to $2 n^2 - 2 \sum_{1 \leq i, j \leq n} X_i X_j$, and we equivalently care about the concentration properties of $\sum_{1 \leq i, j \leq n} X_i X_j$. This is just the simplest possible Rademacher chaos of order two.
Now, you might try looking for lower bounds on the probability that the 2nd order Rademacher chaos is outside some radius $t$ of its expectation. Upper bounds are known, for instance in http://projecteuclid.org/download/pdf_1/euclid.aop/1055425791 and the references in that paper to earlier results by Talagrand and by Ledoux. Perhaps they've also shown optimality of the rates, which would be sufficient for your purposes. I don't think the linked paper shows rate optimality, but try the Talagrand or Ledoux papers.
A: Wouldn't it suffice simply to calculate the dispersion of $f$? You can calculate it easily: each term is equal to $\langle X_i-X_j, X_i-X_j \rangle$, so it's quite easy to work with them, there are no square roots or other nasties. And to calculate the dispersion you need only individual dispersions (easy), and covariances (not too difficult, too). 
(Perhaps, you'll find it easier to use linearity first, so you have the sum of the terms of the types $\langle X_i, X_i \rangle$ and $\langle X_i, X_j \rangle$, this will simplify the further calculations.)
Either the dispersion decreases quicker than $1/n$, and then you are done by Chebyshev inequality. Or it does not, and then no hope for quicker than $O(1/\sqrt{n})$ decrease.
My personal opinion here is that the decrease is indeed $O(1/\sqrt{n})$. For the following reason: the number you evaluate is some kind of a variance, evaluated for your sample $X_1,\dots,X_n$. I would be surprised if there existed a method more efficient than the standard $\hat{\sigma}_n^2$, that has a variance of $\sim 1/n$. In fact, my intuition tells me that your $f$ is proportional to $\hat{\sigma}_n^2$ (or to the sum of its diagonal elements in the $d$-dimensional case).
