# Heterogeneous Distributed Big Data Clustering on Sparse Grids
## 3 Estimating Densities on Sparse Grids
### The Sparse Grid Density Estimation  
The sparse grid density estimation, originally proposed by [PERSON] et al. (2018), uses an initial density guess \(f_{\epsilon}\) than is smoothed using a spline-smoothing approach:  
\[\hat{f}=\operatorname*{arg\,min}_{u\in V}\int_{\Omega}(u(\mathbf{x})-f_{ \epsilon}(\mathbf{x}))^{2}-D\mathbf{x}+\lambda||Lu||_{L_{2}}^{2}. \tag{9}\]  
This approach results in a function \(\hat{f}\in V\) that balances closeness to the initial density guess with the regularization term \(||Lu||_{L_{2}}^{2}\) that enforces smoothness on the resulting density function. The regularization parameter \(\lambda\) controls the degree of smoothness of \(\hat{f}\). \(L\) usually is some differential operator. We use the initial density guess proposed by [PERSON] et al. that places a Dirac delta function \(\delta_{\mathbf{x}_{i}}\) at every data point \(\mathbf{x}_{i}\):  
\[f_{\epsilon}:=\frac{1}{m}\sum_{i=1}^{m}\delta_{\mathbf{x}_{i}}. \tag{10}\]  
As in prior work, we compute the best sparse grid function \(u\in V_{n}^{(1)}\) and use a surplus-based regularization approach [PERSON] et al. (2018). Therefore, the problem to solve is  
\[\hat{f}=\operatorname*{arg\,min}_{u\in V_{n}^{(1)}}\int_{\Omega}(u(\mathbf{x} )-f_{\epsilon}(\mathbf{x}))^{2}-D\mathbf{x}+\lambda\sum_{i=1}^{N}\alpha_{i}^{ 2}. \tag{11}\]  
This formulation leads to a system of linear equations  
\[(B+\lambda I)\boldsymbol{\alpha}=\mathbf{b}, \tag{12}\]  
with \(B_{ij}=(\phi_{i},\phi_{j})_{L_{2}}\), the identity matrix \(I\) and \(b_{i}=\frac{1}{m}\sum_{j=1}^{m}\phi_{i}(\mathbf{x}_{j})\).  
We solve this system of linear equations with a conjugate gradient solver (CG). Given this iterative solver, two major operations need to be performed: calculating the right-hand side once and computing the matrix-vector product \(\mathbf{v}^{\prime}=(B+\lambda I)\mathbf{v}\) in every CG iteration. The calculation of the right-hand side is straightforward. However, the matrix-vector product requires efficient computations of the \(L_{2}\) inner product of pairs of basis functions:  
\[(\phi_{\mathbf{1},\mathbf{i}},\phi_{\mathbf{I}^{\prime},\mathbf{ i}^{\prime}})_{L_{2}} =\int_{\Omega}\phi_{\mathbf{1},\mathbf{i}}(\mathbf{x})\phi_{ \mathbf{I}^{\prime},\mathbf{i}^{\prime}}(\mathbf{x})d\mathbf{x} \tag{13}\] \[=\int_{0}^{1}\phi_{\mathbf{1}_{1},\mathbf{i}_{1}}(x_{1})\phi_{ \mathbf{i}_{1}^{\prime},\mathbf{i}_{1}^{\prime}}(x_{1})dx_{1}\cdot\cdot\cdot \int_{0}^{1}\phi_{\mathbf{i}_{d},\mathbf{i}_{d}}(x_{d})\phi_{\mathbf{i}_{d}^{ \prime},\mathbf{i}_{d}^{\prime}}(x_{d})dx_{d}. \tag{14}\]  
The 1-D integrals can be computed directly:  
\[\int_{0}^{1}\phi_{\mathbf{I},i}(x)\phi_{\mathbf{I}^{\prime},\mathbf{i}^{ \prime}}(x)dx=\begin{cases}\frac{2}{3}h_{\mathbf{I}},&x_{li}=x_{\mathbf{I}^{ \prime},\mathbf{i}^{\prime}},\\ h_{\mathbf{I}^{\prime}}\phi_{\mathbf{I},i}(x_{\mathbf{I}^{\prime}})+h_{ \mathbf{I}}\phi_{\mathbf{I}^{\prime}}(x_{li}),&\text{else}.\end{cases} \tag{15}\]  
We note that in many instances the integral will be zero due to the non-overlapping support of the hat functions.  
Figure 4 shows the effect of varying the regularization parameter \(\lambda\) for a 2-D dataset. \(\lambda\) has to be chosen with care, as too small values might split a single cluster into multiple clusters. On the other hand, if \(\lambda\) is too large separate clusters could be part of the same high-density region.