To understand the intriguing many-body states and effects in the correlated quantum materials, inference of the microscopic effective Hamiltonian from experiments constitutes an important yet very challenging inverse problem. Here we propose an unbiased and efficient approach learning the effective Hamiltonian through the many-body analysis of the measured thermal data. Our approach combines the strategies including the automatic gradient and Bayesian optimization with the thermodynamics many-body solvers including the exact diagonalization and the tensor renormalization group methods. We showcase the accuracy and powerfulness of the Hamiltonian learning by applying it firstly to the thermal data generated from a given spin model, and then to realistic experimental data measured in the spin-chain compound copper nitrate and triangular-lattice magnet TmMgGaO4. The present automatic approach constitutes a unified framework of many-body thermal data analysis in the studies of quantum magnets and strongly correlated materials in general.