The usual formulation of the least squares collocation method requires the inversion of matrices with dimensions equal to the number of observations. Such a procedure poses difficult numerical problems when large numbers of observations are involved. It is shown in this report how this difficulty can be overcome. The size of the matrices to be inverted can be reduced to the number of unknown quantities, i.e. the amount of numerical work will be the same as for the corresponding adjustment problem. Numerical questions as the treatment of orbital parameters, the determination of the speed of convergence, the stability of different solutions, and the choice of a covariance function are discussed. The determination of tesseral harmonic coefficients from a combination of satellite and gravimetric data is treated in some detail and the results are related to previous work.