I do not know if you are still interested in solving this issue. If this is the case, I would suggest to use a different approach if it is ok with you...anyway this answer might be useful to somebody else in the future (I hope at least).
An option to smooth multivariate histograms, is to use P-splines and fit the array of counts as suggested in the comment at your question (see this reference for example). P-splines combine B-spline bases and finite difference penalties (see Eilers and Marx, 1991).
To keep the notation simple, suppose that we want to smooth a 2d hist (the method generalizes to more dim as you will read in the reference I mentioned above). Let define $Y$ the 2-dim array of counts collected on a grid of $n_{1} \times n_{2}$ grid of bins. We can model the expected counts as Poisson using tensor product B-splines as:
$$
\mu = E[\mbox{vec}(Y)] = \exp[(B_{2} \otimes B_{1}) \eta]
$$
where $B_{1}$ and $B_{2}$ are $(n_{1} \times m_{1})$ and $(n_{2} \times m_{2})$ B-spline matrices and $\eta$ is a $m_{1}m_{2}$ vector of unknown spline coefficients to be estimated from the data. The optimal spline coefficients can be estimated by maximizing the following penalized log-likelihood (foe example by using IWLS):
$$
\ell = \sum\log(\mu^{\mbox{vec(Y)}} \exp(-\mu)) - \lambda (\eta^{\top}P\eta)
$$
Again, to keep notation simple, I will consider the isotropic 2d penalty
$$
P = \lambda(I_{n_{2}} \otimes \Delta_{1}^{\top} \Delta_{1} + \Delta_{2}^{\top} \Delta_{2} \otimes I_{n_{1}})
$$
where $\Delta$ is the finite difference matrix operator (I will use here 2nd order differences) and the $\lambda$ parameter regulates the amount of smoothness required for the final estimates (please note that here I consider a single $\lambda$...it is possible to have 2 of them in the definition of $P$ above).
Below you will find a small code reproducing an example (I left some comment in the R-code to help readability, I hope). In this case, I will use a "2D Whittaker smoother" which is a special case of P-spline (you can read about the relationship between these two methods in Eilers (2003) for example).
rm(list = ls())
library(colorout)
library(gplots)
library(fields)
# Simulate some data
set.seed(1234)
N = 500
x = 5 + rnorm(N) * 0.5
y = 2 + rnorm(N) * 0.5
# Bi-variate histogram
h = hist2d(data.frame(x, y), nbins = 25, show = FALSE)
# Data handling
x_mid = h$x
y_mid = h$y
cnt = h$counts
cnt_vec = as.vector(cnt)
n = dim(cnt)[1]
m = dim(cnt)[2]
# Bases stuffs (Whittaker) & penalties (isotropic)
dd = 2
lm = 1e2 # Fixed lambda - can be selected
E1 = diag(n)
E2 = diag(m)
Eo = E1 %x% E2
P1 = crossprod(diff(E1, diff = dd))
P2 = crossprod(diff(E2, diff = dd))
Po = lm * diag(m) %x% P1 + lm * P2 %x% diag(n) # isotropic smoothing, can be made anisotropic
# Smoothing
tol = 1e-8
eta = log(cnt_vec + 1)
for (it in 1:50)
{
mu = exp(eta)
z = cnt_vec - mu + mu * eta
a = solve(t(Eo) %*% (c(mu) * Eo) + Po, t(Eo) %*% z)
etnew = Eo %*% a
de = max(abs(etnew - eta))
cat('Crit', it, de, '\n')
if(de < tol) break
eta = etnew
}
cnt_fit = matrix(exp(eta), nc = m, nr = n)
# Plot results
image.plot(x_mid, y_mid, cnt)
contour(x_mid, y_mid, cnt_fit, add = TRUE, lwd = 2, nlevels = 25)
Finally, I hope I got the notation right and that my answer replies (at least partially) the original question.
