A data-set of nearly 100,000 symmetry unique multi-configurational ab initio points for methane were generated at the (AE)-MRCI-F12(Q)/CVQZ-F12 level, including energies beyond 30,000 cm−1 above the minimum and fit into potential energy surfaces (PESs) by several permutation invariant schemes. A multi-expansion interpolative fit combining interpolating moving least squares (IMLS) fitting and permutation invariant polynomials (PIP) was able to fit the complete data-set to a root-mean-square deviation of 1.0 cm−1 and thus was used to benchmark the other fitting methods. The other fitting methods include a single PIP expansion and two neural network (NN) based approaches, one of which combines NN with PIP. Full-dimensional variational vibrational calculations using a contracted-iterative method (and a Lanczos eigensolver) were used to assess the spectroscopic accuracy of the electronic structure method. The results show that the NN-based fitting approaches are able to fit the data-set remarkably accurately with the PIP-NN method producing levels in remarkably close agreement with the PIP-IMLS benchmark. The (AE)-MRCI-F12(Q)/CVQZ-F12 electronic structure method produces vibrational levels of near spectroscopic accuracy and a superb equilibrium geometry. The levels are systematically slightly too high, beginning at ∼ 1-2 cm−1 above the fundamentals and becoming correspondingly higher for overtones. The PES is therefore suitable for small ab initio or empirical corrections and since it is based on a multi-reference method, can be extended to represent dynamically relevant dissociation channels.