kk and ∆ (l) (k-1)k are given in the Supporting Information. Being the algorithm stochastic in its nature, it is of interest to analyze the distribution of the energy estimations in terms of the number of random samples (R) as well as the dimension of the Krylov subspace identifying the number of Lanczos quadrature points. This is shown in Figure where three distributions relative to three different values of R are reported for a Krylov subspace of dimension 20. The SL distributions are centered close to zero and this denotes a negligible bias toward the target energy value obtaineded via a direct diagonalization. Moreover, the distributions narrow with the increment of R as one would expect, thus allowing for a systematic improvement of the estimation. We note that by decreasing the dimension of the employed Krylov subspace to 10, the distributions remain well centered around zero with a bias smaller than 1% on the total MBD energy even for low values of R.