Computation of covariance matrices from observed data is an important problem, as such matrices are used in applications such as PCA, LDA, and increasingly in the learning and application of probabilistic graphical models. However, computing an empirical covariance matrix is not always an easy problem.There are two key difficulties associated with computing such a matrix from a very high-dimensional data set. The first problem is over-fitting. For a p-dimensional covariance matrix, there are p(p − 1)/2 unique, off-diagonal entries in the empirical covariance matrixŜ; for large p (say, p > 10 5 ) the size n of the data set is often much smaller than the number of covariances to compute. Over-fitting is a concern in any situation where the number of parameters learned can greatly exceed the size of the data set. Thus, there are strong theoretical reasons to expect that for highdimensional data-even Gaussian data-the empirical covariance matrix is not a good estimate for the true covariance matrix underlying the generative process.The second problem is computational. Computing a covariance matrix takes O(np 2 ) time. For large p (greater than 10,000) and n much greater than p, this is debilitating.The first problem (over-fitting) has been studied in depth in both the statistics and machine learning literature, but the second problem (computation) has been ignored, presumably because these fields have typically been concerned with relatively small data sets. In this thesis, we consider how both of these difficulties can be handled simultaneously. Specifically, a key regularization technique for high-dimensional covariance estimation is thresholding, where the smallest or least significant entries in the covariance matrix are simply dropped and replaced with the value 0. This suggests an obvious way to address the computational difficulty as well: first, compute the identities of the K entires in the covariance matrix that are actually important in the sense that they will not be removed during thresholding, and then in a second step, compute the values of those entries. This can be done in O(Kn) time. If K << p 2 and the identities of the important entries can be computed in reasonable time, then this is a big win.The key technical contribution of this thesis is the design and implementation of two different distributed algorithms for approximating the identities of the important entries quickly, using sampling. We have implemented these methods and tested them using an 800 core compute cluster. Experiments have been run using real data sets having millions of data points and up to 40, 000 dimensions. These experiments show that the proposed methods are both accurate and efficient.