  The non-parametric estimation of covariance lies at the heart of functional
data analysis, whether for curve or surface-valued data. The case of a
two-dimensional domain poses both statistical and computational challenges,
which are typically alleviated by assuming separability. However, separability
is often questionable, sometimes even demonstrably inadequate. We propose a
framework for the analysis of covariance operators of random surfaces that
generalises separability, while retaining its major advantages. Our approach is
based on the expansion of the covariance into a series of separable terms. The
expansion is valid for any covariance over a two-dimensional domain. Leveraging
the key notion of the partial inner product, we extend the power iteration
method to general Hilbert spaces and show how the aforementioned expansion can
be efficiently constructed in practice. Truncation of the expansion and
retention of the leading terms automatically induces a non-parametric estimator
of the covariance, whose parsimony is dictated by the truncation level. The
resulting estimator can be calculated, stored and manipulated with little
computational overhead relative to separability. Consistency and rates of
convergence are derived under mild regularity assumptions, illustrating the
trade-off between bias and variance regulated by the truncation level. The
merits and practical performance of the proposed methodology are demonstrated
in a comprehensive simulation study and on classification of EEG signals.
