The Nystrom Method is a sampling-based or CUR-like low-rank approximation to symmetric positive semi-definite (SPSD) matrices such as Laplacians and kernels that arise in machine learning. The method has received a great deal of attention in recent years, much of which has focused on the use of randomized algorithms to compute efficiently the low-rank approximation. Motivated by conflicting and contradictory claims in the literature, we have performed a detailed empirical evaluation of the performance quality and running time of sampling-based and projection-basedlow-rank approximation methods on a diverse suite of SPSD matrices. Our main conclusions are threefold: first, our results highlight complementary aspects of projection methods versus sampling based on the statistical leverage scores; second, our results elucidate the effect of popular "design decisions" on the structural properties of the matrices and thus on the applicability of different approximation algorithms; and third, our results show that sampling-based algorithms that use provably-accurateapproximate leverage scores can have running times comparable to random projection algorithms. In addition, our empirical results demonstrate that prior theory was so weak as to not even be a qualitative guide to practice. Thus, we complement our empirical results with a suite of worst-case theoretical bounds for both random sampling and random projections methods. These bounds are qualitatively superior to existing bounds, e.g., improved additive-error bounds for the spectral and Frobenius norm errors and relative-error bounds for the trace norm error, and they can be used to understand the complementary aspects of projection versus sampling-based methods. We will also describe how recent numerical implementations of random sampling and random projection algorithms can be used in this setting in much larger-scale applications.