Proof that Generators of Centroidal Voronoi Tesselation are Given by Conditional Mean I am trying to understand a proof in "Centroidal voronoi
tessellations: applications and algorithms" from Du et al. In particular, I do not understand Proposition 3.1, stating that a necessary condition for an "optimal" Voronoi tesselation its generators need be equal to the centroids of the cell they define.
It begins with a first order variation of the function F, then we should divide by $\epsilon$ and construct the limes. Unfortunately, I do not understand how we get this way from the formula for the first order variation to $z_j$.
Copied from the text, the step I do not understand is going from:
$F(z_j + \epsilon v) - F(z_j) = \int_{y \in Y_j} p(y) |y - z_j - \epsilon v|^2 - |y - z_j|^2 dy$
to 
$z_j = \frac{\int_{y \in Y_j} yp(y)dy}{\int_{y \in Y_j} p(y) dy}$.
 A: Expanding the norm squared yields $\int_{y\in V_j} \rho({\bf y})\left( \left|{\bf y} -{\bf z}_j - \epsilon {\bf v}\right|^2 -  \left|{\bf y} -{\bf z}_j\right|^2 \right) d{\bf y}$
\begin{align}
& = \int_{y\in V_j} \rho({\bf y})\left( \left|{\bf y-{\bf z}_j}\right|^2 + \epsilon^2 \left|{\bf v}\right|^2 -2\epsilon \left({\bf y-{\bf z}_j}\right)\cdot{\bf v} - \left|{\bf y-{\bf z}_j}\right|^2\right) d{\bf y} \\
& = \int_{y\in V_j} \rho({\bf y}) \left( \epsilon^2 \left|{\bf v}\right|^2 -2\epsilon \left({\bf y-{\bf z}_j}\right)\cdot{\bf v}\right) d{\bf y}.
\end{align}
When dividing by $\epsilon$ and taking the limit, the $\epsilon^2$ term drops out giving,
\begin{align}
\lim_{\epsilon \rightarrow 0} \frac{{\mathcal F}({\bf z}_j + \epsilon {\bf v}) - {\mathcal F}({\bf z}_j) }{\epsilon} & = -2 \int_{y\in V_j} \rho({\bf y}) \left({\bf y-{\bf z}_j}\right)\cdot{\bf v} d{\bf y}\\
 & = -2{\bf v} \cdot \int_{y\in V_j} \rho({\bf y}) \left({\bf y}-{\bf z}_j\right) d{\bf y}.
\end{align}
If ${\bf z}_j$ are minimizers of this functional, then this derivative must be zero for any possible vector ${\bf v}$. Thus 
$$
{\bf 0} = -2 \int_{y\in V_j} \rho({\bf y}) \left({\bf y-{\bf z}_j}\right) d{\bf y}.
$$
Rearranging terms gives,
\begin{align}
\int_{y\in V_j} \rho({\bf y}){\bf y} d{\bf y} & = \int_{y\in V_j} \rho({\bf y}){\bf z}_j d{\bf y} \\
& = \left(\int_{y\in V_j} \rho({\bf y}) d{\bf y} \right){\bf z}_j.
\end{align}
Dividing by $\int_{y\in V_j} \rho({\bf y}) d{\bf y}$ gives the expression reported in the paper,
$$
{\bf z}_j = \frac{ \int_{y\in V_j} \rho({\bf y}){\bf y} d{\bf y} } { \int_{y\in V_j} \rho({\bf y}) d{\bf y} }
$$
