We present a computational approach to high order matching of datasets in \mathbb{R}^{d}. Those are matchings based on data affinity measures that score the matching of more than two pairs of points at a time. High order affinities are represented by tensors and the matching is then given by a rank-one approximation of the affinity tensor and a corresponding discretization. Our approach is rigorously justified by extending Zass and Shashua's hypergraph matching [40] to high order spectral matching. This paves the way for a computationally efficient dual marginalization-spectral matching scheme. We also show that based on the spectral properties of random matrices, affinity tensors can be randomly sparsified while retaining the matching accuracy. Our contributions are experimentally validated by applying them to synthetic as well as real datasets.