  The derivation of cosmological parameters from astrophysical data sets
routinely involves operations counts which scale as O(N^3) where N is the
number of data points. Currently planned missions, including MAP and Planck,
will generate sky maps with N_d = 10^6 or more pixels. Simple ``brute force''
analysis, applied to such mega-pixel data, would require years of computing
even on the fastest computers. We describe an algorithm which allows estimation
of the likelihood function in the direct pixel basis. The algorithm uses a
conjugate gradient approach to evaluate chi-squared and a geometric
approximation to evaluate the determinant. Monte Carlo simulations provide a
correction to the determinant, yielding an unbiased estimate of the likelihood
surface in an arbitrary region surrounding the likelihood peak. The algorithm
requires O(N_d^{3/2}) operations and O(N_d) storage for each likelihood
evaluation, and allows for significant parallel computation.
