In one-dimensional quantum systems with short-range interactions, a set of leading numerical methods is based on matrix product states, whose bond dimension determines the amount of computational resources required by these methods. We prove that a thermal state at constant inverse temperature $beta$ has a matrix product representation with bond dimension $e^{tilde O(sqrt{betalog(1/epsilon)})}$ such that all local properties are approximated to accuracy $epsilon$. This justifies the common practice of using a constant bond dimension in the numerical simulation of thermal properties.