Write $X_{B,b} = \{\alpha \in \mathbb{Z}^B : \sum_j \alpha_j = b\}$. Now using the convention $0^0 \equiv 1$, define the matrix $W_{\alpha, \alpha'} := \alpha^{\alpha'}$. For arbitrary $f:X_{B,b} \rightarrow \mathbb{K}$ we can write $f_\alpha \equiv \sum_{\alpha'} c_{\alpha'}$ from which it follows that $c = W^{-1}f$. The facts that this procedure is well-defined, and that $W$ possesses an inverse, follow from a result in multivariate interpolation assuring us that the Lagrange interpolation problem on $X_{B,b}$ is "poised".
Actually, although generic discrete point sets admit a specific multivariate Lagrange interpolation protocol that satisfies many desirable properties, only $X_{B,b}$ does it so beautifully. As a result, we obtain a Lagrange interpolation: $f_{\mathfrak{I}}(x) := \sum_\alpha (W^{-1}f)_\alpha x^\alpha$ which satisfies $f_{\mathfrak{I}}(\alpha) = f(\alpha)$.
You can use this to define differencing schemes on the hexagonal grid by considering $B = 3$. An example of the interpolation is shown.
Define $d_{\mathfrak{I}} f := d(f_{\mathfrak{I}})|_{X_{B,b}}$. Note (e.g.) that
$\partial_j f_{\mathfrak{I}} = \sum_\alpha (W^{-1} f)_\alpha \partial_j x^\alpha = \sum_\alpha (W^{-1} f)_\alpha \frac{\alpha_j}{x_j} \partial_j x^\alpha$
(suitably interpreted) is easy to compute in silico. Explicitly, set
$\left(W_{(\partial_j)}\right)_{\alpha, \alpha'} := \frac{\alpha'_j}{\alpha_j} \alpha^{\alpha'}, \quad \left(\mathcal{W}_{(\partial_j)}\right)_{x, \alpha'} := \frac{\alpha'_j}{x_j} x^{\alpha'}.$
Then
$\partial_j f_{\mathfrak{I}} = \mathcal{W}_{(\partial_j)} W^{-1} f, \quad \partial_j f \equiv W_{(\partial_j)} W^{-1} f.$