  We propose a fast penalized spline method for bivariate smoothing. Univariate
P-spline smoothers (Eilers and Marx, 1996) are applied simultaneously along
both coordinates. The new smoother has a sandwich form which suggested the name
"sandwich smoother" to a referee. The sandwich smoother has a tensor product
structure that simplifies an asymptotic analysis and it can be fast computed.
We derive a local central limit theorem for the sandwich smoother, with simple
expressions for the asymptotic bias and variance, by showing that the sandwich
smoother is asymptotically equivalent to a bivariate kernel regression
estimator with a product kernel. As far as we are aware, this is the first
central limit theorem for a bivariate spline estimator of any type. Our
simulation study shows that the sandwich smoother is orders of magnitude faster
to compute than other bivariate spline smoothers, even when the latter are
computed using a fast GLAM (Generalized Linear Array Model) algorithm, and
comparable to them in terms of mean squared integrated errors. We extend the
sandwich smoother to array data of higher dimensions, where a GLAM algorithm
improves the computational speed of the sandwich smoother. One important
application of the sandwich smoother is to estimate covariance functions in
functional data analysis. In this application, our numerical results show that
the sandwich smoother is orders of magnitude faster than local linear
regression. The speed of the sandwich formula is important because functional
data sets are becoming quite large.
