  We construct interpolation operators 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.
