Analytic representations of coupled potential energy surfaces (PESs) of polyatomic molecules are valuable resources because one can perform detailed dynamics simulation of the corresponding system with appropriately large ensembles, ultralong simulation times, accurate electronic structure (which would be unaffordable without fitting), and accurate dynamics algorithms. The development of PES-learning methodologies has enabled highly accurate fitting of coupled global PESs for complex systems and systems with dense manifolds of states. Critical aspects of an accurate learning of a set of global coupled PESs include: (1) The fit should maintain -exactly or to a practically useful approximation -the invariance properties of potential energy surfaces, namely the invariance with respect to the translation, rotation, and permutation of identical nuclei. (2) The PESs should become flat in asymptotic regions, i.e., when a subsystem is dissociated from the rest of the system. (3) For coupled PESs in the regions where adiabatic PESs cross, i.e. conical intersections, the learned surfaces should provide the correct topology of the conical intersections without the potential error of double crossing.