Hierarchical equations of motion theory for Drude dissipation is optimized, with a convenient convergence criterion proposed in advance of numerical propagations. The theoretical construction is on the basis of a Padé spectrum decomposition that has been qualified to be the best sum-over-poles scheme for quantum distribution function. The resulting hierarchical dynamics under the a priori convergence criterion are exemplified with a benchmark spin-boson system, and also the transient absorption and related coherent two-dimensional spectroscopy of a model exciton dimer system. We combine the present theory with several advanced techniques such as the block hierarchical dynamics in mixed Heisenberg-Schrödinger picture and the on-the-fly filtering algorithm for the efficient evaluation of third-order optical response functions.