Many gauge field theories can be described using a multisymplectic Lagrangian formulation, where the configuration manifold is the space of Lorentzian metrics. Group-equivariant interpolation spaces are critical to the construction of geometric structure-preserving discretizations of such problems, since they can be used to construct a variational discretization that exhibits a discrete Noether's theorem. We approach this problem more generally, by considering interpolation spaces for functions taking values in a symmetric space --- a smooth manifold with an inversion symmetry about every point.

Key to our construction is the observation that every symmetric space can be realized as a homogeneous space whose cosets have canonical representatives by virtue of the generalized polar decomposition -- a generalization of the well-known factorization of a real nonsingular matrix into the product of a symmetric positive-definite matrix times an orthogonal matrix. By interpolating these canonical coset representatives, we derive a family of structure-preserving interpolation operators for symmetric space-valued functions. As applications, we construct interpolation operators for the space of Lorentzian metrics, the space of symmetric positive-definite matrices, and the Grassmannian. In the case of Lorentzian metrics, our interpolation operators provide a family of finite elements for numerical relativity that are frame-invariant and have signature which is guaranteed to be Lorentzian pointwise. We illustrate their potential utility by interpolating the Schwarzschild metric numerically.