We develop the multi-layer multi-configuration time-dependent Hartree method for bosons (ML-MCTDHB), a variational numerically exactab initiomethod for studying the quantum dynamicsand stationary properties of general bosonic systems. ML-MCTDHB takes advantage of the per-mutation symmetry of identical bosons, which allows for investigations of the quantum dynamicsfrom few to many-body systems. Moreover, the multi-layer feature enables ML-MCTDHB to de-scribe mixed bosonic systems consisting of arbitrary many species. Multi-dimensional as well asmixed-dimensional systems can be accurately and efficiently simulated via the multi-layer expan-sion scheme. We provide a detailed account of the underlying theory and the corresponding imple-mentation. We also demonstrate the superior performance by applying the method to the tunnelingdynamics of bosonic ensembles in a one-dimensional double well potential, where a single-speciesbosonic ensemble of various correlation strengths and a weakly interacting two-species bosonic en-semble are considered.