||Quantum dissipation theory (QDT) plays a key role describing the reduced system dynamics under influence of its environment. Conventional QDTs, i.e, quantum master equations (QMEs), explore weak system-bath coupling and also often Markovian approximations. However, there are many of problem in mesoscopic systems, where perturbative and Markovian theories would not be adequate and lead to even physical results. The existing exact theories such as Feynman-Vernon influence functional path integral formalism are numerically too expensive to realistic systems. In this thesis, I focus on the hierarchical equations of motion (HEOM) approach, which is mathematically equivalent to Feynman-Vernon influence functional theory, but numerically and operationally much more implementable. The dynamical variables in HEOM are the reduced system operator and its associated auxiliary density operators (ADOs). The explicit hierarchical construction depends on the way of decomposing the interacting environment correlation functions into their memory frequency components. This is the statistical environment basis set problem in HEOM. On the other hand, the converged level of hierarchy is dictated physically by the combined effects of memory, system-bath coupling strength, and anharmonicity. The total number of ADOs, as it depends on the environment basis set size and the level of hierarchy, resembles the full configuration interaction problem in many-particle quantum mechanism. The achievement of this thesis research is mainly related to the three important issues in the choice of statistical environment basis set for efficient HEOM construction. The first one is to identify the best basis set for spanning over the stochastic bath space. The second one concerns the minimum basis set. It requests to have not only the smallest number of decomposition terms, but also a priori accuracy control criterion on the resulting HEOM dynamics for any given system under bath influence at finite temperature. The third issue is about the possibility of at least partial inclusion of the off-basis-set residue effect on the HEOM dynamics. These issues are all related to the memory frequency decomposition of environment correlation functions that satisfy the fermionic or bosonic fluctuation-dissipation theory (FDT). To enable numerical study on realistic systems, which are often non-Markovian or non-perturbative, an optimized HEOM with minimum but physically sufficient basis sets is crucial. To achieve that, a family of three Padé spectrum decomposition (PSD) schemes are proposed to replace the existing Mastubara spectrum decomposition (MSD), which minimizes the number of decomposition terms. Besides, a very efficient on-the-fly filter algorithm is also applied for the numerical HEOM propagation. Finally, serval of numerical applications are presented for model systems study. These include thermal driven charge transport, high-order quantum transport, and time-dependent transient current.