  The paper motivates high dimensional smoothing with penalized splines and its
numerical calculation in an efficient way. If smoothing is carried out over
three or more covariates the classical tensor product spline bases explode in
their dimension bringing the estimation to its numerical limits. A recent
approach by Siebenborn and Wagner(2019) circumvents storage expensive
implementations by proposing matrix-free calculations which allows to smooth
over several covariates. We extend their approach here by linking penalized
smoothing and its Bayesian formulation as mixed model which provides a
matrix-free calculation of the smoothing parameter to avoid the use of
high-computational cross validation. Further, we show how to extend the ideas
towards generalized regression models. The extended approach is applied to
remote sensing satellite data in combination with spatial smoothing.
