Generalized least square problems with non-diagonal weights arise frequently in an estimation of two dimensional images from data of cosmological as well as astro- or geophysical observations. As the observational data sets keep growing at Moore's rate, with their volumes exceeding tens and hundreds billions of samples, the need for fast and efficiently parallelizable iterative solvers is generally recognized. In this work we study performance of two-level preconditioners in the context of iterative solvers of the generalized least square systems, where the weights are assumed to be described by a block-diagonal matrix with Toeplitz blocks. Such cases are physically well motivated and arise whenever the instrumental noise displays a piece-wise stationary behavior. Our iterative algorithm is based on a conjugate gradient method with a parallel two-level preconditioner (2lvl-PCG) for which we construct its coarse space from a limited number of sparse vectors estimated solely from coefficients of the initial linear system. Our prototypical application is the map-making problem in the Cosmic Microwave Background observations. We show experimentally that our parallel implementation of 2lvl-PCG outperforms by as much as a factor 5 standard one-level PCG in terms of both the convergence rate and time to the solution.

Authors: M. Szydlarski, L. Grigori and R. Stompor