Scalable iterative methods for sampling from massive Gaussian random vectors

Sampling from Gaussian Markov random fields (GMRFs), that is multivariate Gaussian ran- dom vectors that are parameterised by the inverse of their covariance matrix, is a fundamental problem in computational statistics. In this paper, we show how we can exploit arbitrarily accu- rate approximations to a GMRF to speed up Krylov subspace sampling methods. We also show that these methods can be used when computing the normalising constant of a large multivariate Gaussian distribution, which is needed for both any likelihood-based inference method. The method we derive is also applicable to other structured Gaussian random vectors and, in particu- lar, we show that when the precision matrix is a perturbation of a (block) circulant matrix, it is still possible to derive O(n log n) sampling schemes.Comment: 17 Pages, 4 Figure

Accelerating Cosmic Microwave Background map-making procedure through preconditioning

Estimation of the sky signal from sequences of time ordered data is one of the key steps in Cosmic Microwave Background (CMB) data analysis, commonly referred to as the map-making problem. Some of the most popular and general methods proposed for this problem involve solving generalised least squares (GLS) equations with non-diagonal noise weights given by a block-diagonal matrix with Toeplitz blocks. In this work we study new map-making solvers potentially suitable for applications to the largest anticipated data sets. They are based on iterative conjugate gradient (CG) approaches enhanced with novel, parallel, two-level preconditioners. We apply the proposed solvers to examples of simulated non-polarised and polarised CMB observations, and a set of idealised scanning strategies with sky coverage ranging from nearly a full sky down to small sky patches. We discuss in detail their implementation for massively parallel computational platforms and their performance for a broad range of parameters characterising the simulated data sets. We find that our best new solver can outperform carefully-optimised standard solvers used today by a factor of as much as 5 in terms of the convergence rate and a factor of up to $4$ in terms of the time to solution, and to do so without significantly increasing the memory consumption and the volume of inter-processor communication. The performance of the new algorithms is also found to be more stable and robust, and less dependent on specific characteristics of the analysed data set. We therefore conclude that the proposed approaches are well suited to address successfully challenges posed by new and forthcoming CMB data sets.Comment: 19 pages // Final version submitted to A&

Accelerating Cosmic Microwave Background map-making procedure through preconditioning

Estimation of the sky signal from sequences of time ordered data is one of the key steps in Cosmic Microwave Background (CMB) data analysis, commonly referred to as the map-making problem. Some of the most popular and general methods proposed for this problem involve solving generalised least squares (GLS) equations with non-diagonal noise weights given by a block-diagonal matrix with Toeplitz blocks. In this work we study new map-making solvers potentially suitable for applications to the largest anticipated data sets. They are based on iterative conjugate gradient (CG) approaches enhanced with novel, parallel, two-level preconditioners. We apply the proposed solvers to examples of simulated non-polarised and polarised CMB observations, and a set of idealised scanning strategies with sky coverage ranging from nearly a full sky down to small sky patches. We discuss in detail their implementation for massively parallel computational platforms and their performance for a broad range of parameters characterising the simulated data sets. We find that our best new solver can outperform carefully-optimised standard solvers used today by a factor of as much as 5 in terms of the convergence rate and a factor of up to $4$ in terms of the time to solution, and to do so without significantly increasing the memory consumption and the volume of inter-processor communication. The performance of the new algorithms is also found to be more stable and robust, and less dependent on specific characteristics of the analysed data set. We therefore conclude that the proposed approaches are well suited to address successfully challenges posed by new and forthcoming CMB data sets.Comment: 19 pages // Final version submitted to A&

Spectral behavior of preconditioned non-Hermitian multilevel block Toeplitz matrices with matrix-valued symbol

This note is devoted to preconditioning strategies for non-Hermitian multilevel block Toeplitz linear systems associated with a multivariate Lebesgue integrable matrix-valued symbol. In particular, we consider special preconditioned matrices, where the preconditioner has a band multilevel block Toeplitz structure, and we complement known results on the localization of the spectrum with global distribution results for the eigenvalues of the preconditioned matrices. In this respect, our main result is as follows. Let $I_k:=(-\pi,\pi)^k$, let $\mathcal M_s$ be the linear space of complex $s\times s$ matrices, and let $f,g:I_k\to\mathcal M_s$ be functions whose components $f_{ij},\,g_{ij}:I_k\to\mathbb C,\ i,j=1,\ldots,s,$ belong to $L^\infty$. Consider the matrices $T_n^{-1}(g)T_n(f)$, where $n:=(n_1,\ldots,n_k)$ varies in $\mathbb N^k$ and $T_n(f),T_n(g)$ are the multilevel block Toeplitz matrices of size $n_1\cdots n_ks$ generated by $f,g$. Then $\{T_n^{-1}(g)T_n(f)\}_{n\in\mathbb N^k}\sim_\lambda g^{-1}f$, i.e. the family of matrices $\{T_n^{-1}(g)T_n(f)\}_{n\in\mathbb N^k}$ has a global (asymptotic) spectral distribution described by the function $g^{-1}f$, provided $g$ possesses certain properties (which ensure in particular the invertibility of $T_n^{-1}(g)$ for all $n$) and the following topological conditions are met: the essential range of $g^{-1}f$, defined as the union of the essential ranges of the eigenvalue functions $\lambda_j(g^{-1}f),\ j=1,\ldots,s$, does not disconnect the complex plane and has empty interior. This result generalizes the one obtained by Donatelli, Neytcheva, Serra-Capizzano in a previous work, concerning the non-preconditioned case $g=1$. The last part of this note is devoted to numerical experiments, which confirm the theoretical analysis and suggest the choice of optimal GMRES preconditioning techniques to be used for the considered linear systems.Comment: 18 pages, 26 figure

Preconditioners for Krylov subspace methods: An overview

When simulating a mechanism from science or engineering, or an industrial process, one is frequently required to construct a mathematical model, and then resolve this model numerically. If accurate numerical solutions are necessary or desirable, this can involve solving large-scale systems of equations. One major class of solution methods is that of preconditioned iterative methods, involving preconditioners which are computationally cheap to apply while also capturing information contained in the linear system. In this article, we give a short survey of the field of preconditioning. We introduce a range of preconditioners for partial differential equations, followed by optimization problems, before discussing preconditioners constructed with less standard objectives in mind

FAST SOLUTION METHODS FOR CONVEX QUADRATIC OPTIMIZATION OF FRACTIONAL DIFFERENTIAL EQUATIONS

In this paper, we present numerical methods suitable for solving convex quadratic Fractional Differential Equation (FDE) constrained optimization problems, with box constraints on the state and/or control variables. We develop an Alternating Direction Method of Multipliers (ADMM) framework, which uses preconditioned Krylov subspace solvers for the resulting sub-problems. The latter allows us to tackle a range of Partial Differential Equation (PDE) optimization problems with box constraints, posed on space-time domains, that were previously out of the reach of state-of-the-art preconditioners. In particular, by making use of the powerful Generalized Locally Toeplitz (GLT) sequences theory, we show that any existing GLT structure present in the problem matrices is preserved by ADMM, and we propose some preconditioning methodologies that could be used within the solver, to demonstrate the generality of the approach. Focusing on convex quadratic programs with time-dependent 2-dimensional FDE constraints, we derive multilevel circulant preconditioners, which may be embedded within Krylov subspace methods, for solving the ADMM sub-problems. Discretized versions of FDEs involve large dense linear systems. In order to overcome this difficulty, we design a recursive linear algebra, which is based on the Fast Fourier Transform (FFT). We manage to keep the storage requirements linear, with respect to the grid size $N$, while ensuring an order $N \log N$ computational complexity per iteration of the Krylov solver. We implement the proposed method, and demonstrate its scalability, generality, and efficiency, through a series of experiments over different setups of the FDE optimization problem