The overall computational cost of the SL algorithm is to be mainly addressed to the construction of the tridiagonal matrices ∆ (l) via the Lanczos algorithm which, however, only relies on matrix-vector products and can thus be performed in a linear-scaling fashion. On the other hand, the computational cost required by the diagonalization of ∆ (l) is negligible as the (M + 1) number of quadrature points (and therefore the dimension of the tridiagonal matrix to be diagonalized) is very small, typically not exceeding 50, regardless of the system size. Moreover, the exact diagonalization of small tridiagonal matrices can be performed efficiently by means of several O(N) efficient libraries. We finally note that the SL algorithm is embarrassingly parallelizable with respect to the R random vectors as they involve independent processes thus requiring minimal communication.