Understanding and optimizing the relaxation dynamics of many-body systems is essential both for foundational studies in quantum thermodynamics and for applications such as quantum simulation and quantum computing. Efficient preparation of thermal states of a many-body Hamiltonian is governed by the spectral properties of the associated Lindbladian, in particular its spectral gap, which determines the minimum relaxation rate. In this work, we develop a systematic framework for constructing Lindbladians that prepare thermal states. Our approach reveals a simple relation between the relaxation dynamics at finite and infinite temperatures. The framework is scalable to larger system sizes when implemented using tensor-network methods. We find that efficient thermalization requires that the relaxation dynamics respect the symmetries of the thermal state, which reduces the number of free parameters. By applying gradient-based optimization to the Lindbladians, we enhance the spectral gap and thereby boost thermalization. When applied to both classical and quantum spin models, our method demonstrates a substantial enhancement of the spectral gap. For larger system sizes, our approach provides a variational upper bound and enables a certified lower bound on the spectral gap.