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

Learning detectors quickly using structured covariance matrices

Computer vision is increasingly becoming interested in the rapid estimation of object detectors. Canonical hard negative mining strategies are slow as they require multiple passes of the large negative training set. Recent work has demonstrated that if the distribution of negative examples is assumed to be stationary, then Linear Discriminant Analysis (LDA) can learn comparable detectors without ever revisiting the negative set. Even with this insight, however, the time to learn a single object detector can still be on the order of tens of seconds on a modern desktop computer. This paper proposes to leverage the resulting structured covariance matrix to obtain detectors with identical performance in orders of magnitude less time and memory. We elucidate an important connection to the correlation filter literature, demonstrating that these can also be trained without ever revisiting the negative set

Fast Computation of Smith Forms of Sparse Matrices Over Local Rings

We present algorithms to compute the Smith Normal Form of matrices over two families of local rings. The algorithms use the \emph{black-box} model which is suitable for sparse and structured matrices. The algorithms depend on a number of tools, such as matrix rank computation over finite fields, for which the best-known time- and memory-efficient algorithms are probabilistic. For an \nxn matrix $A$ over the ring \Fzfe, where $f^e$ is a power of an irreducible polynomial f \in \Fz of degree $d$, our algorithm requires \bigO(\eta de^2n) operations in \F, where our black-box is assumed to require \bigO(\eta) operations in \F to compute a matrix-vector product by a vector over \Fzfe (and $\eta$ is assumed greater than \Pden). The algorithm only requires additional storage for \bigO(\Pden) elements of \F. In particular, if \eta=\softO(\Pden), then our algorithm requires only \softO(n^2d^2e^3) operations in \F, which is an improvement on known dense methods for small $d$ and $e$. For the ring \ZZ/p^e\ZZ, where $p$ is a prime, we give an algorithm which is time- and memory-efficient when the number of nontrivial invariant factors is small. We describe a method for dimension reduction while preserving the invariant factors. The time complexity is essentially linear in $\mu n r e \log p,$ where $\mu$ is the number of operations in \ZZ/p\ZZ to evaluate the black-box (assumed greater than $n$) and $r$ is the total number of non-zero invariant factors. To avoid the practical cost of conditioning, we give a Monte Carlo certificate, which at low cost, provides either a high probability of success or a proof of failure. The quest for a time- and memory-efficient solution without restrictions on the number of nontrivial invariant factors remains open. We offer a conjecture which may contribute toward that end.Comment: Preliminary version to appear at ISSAC 201

Recursive-Based PCG Methods for Toeplitz Systems with Nonnegative Generating Functions

In this paper, we consider the solutions of symmetric positive definite, but ill-conditioned, Toeplitz systems An x = b. Here we propose to solve the system by the recursive-based preconditioned conjugate gradient method. The idea is to use the inverse of Am (the principal submatrix of An with the Gohberg--Semencul formula as a preconditioner for An. The inverse of Am can be generated recursively by using the formula until m is small enough. The construction of the preconditioners requires only the entries of An and does not require the explicit knowledge of the generating function f of An. We show that if f is a nonnegative, bounded, and piecewise continuous even function with a finite number of zeros of even order, the spectra of the preconditioned matrices are uniformly bounded except for a fixed number of outliers. Hence the conjugate gradient method, when applied to solving the preconditioned system, converges very quickly. Numerical results are included to illustrate the effectiveness of our approach.published_or_final_versio

A preconditioned MINRES method for nonsymmetric Toeplitz matrices

Circulant preconditioning for symmetric Toeplitz linear systems is well established; theoretical guarantees of fast convergence for the conjugate gradient method are descriptive of the convergence seen in computations. This has led to robust and highly efficient solvers based on use of the fast Fourier transform exactly as originally envisaged in [G. Strang, Stud. Appl. Math., 74 (1986), pp. 171--176]. For nonsymmetric systems, the lack of generally descriptive convergence theory for most iterative methods of Krylov type has provided a barrier to such a comprehensive guarantee, though several methods have been proposed and some analysis of performance with the normal equations is available. In this paper, by the simple device of reordering, we rigorously establish a circulant preconditioned short recurrence Krylov subspace iterative method of minimum residual type for nonsymmetric (and possibly highly nonnormal) Toeplitz systems. Convergence estimates similar to those in the symmetric case are established

Spectral features of matrix-sequences, GLT, symbol, and application in preconditioning Krylov methods, image deblurring, and multigrid algorithms.

The final purpose of any scientific discipline can be regarded as the solution of real-world problems. With this aim, a mathematical modeling of the considered phenomenon is often compulsory. Closed-form solutions of the arising functional equations are usually not available and numerical discretization techniques are required. In this setting, the discretization of an infinite-dimensional linear equation via some linear approximation method, leads to a sequence of linear systems of increasing dimension whose coefficient matrices could inherit a structure from the continuous problem. For instance, the numerical approximation by local methods of constant or nonconstant coefficients systems of Partial Differential Equations (PDEs) over multidimensional domains, gives rise to multilevel block Toeplitz or to Generalized Locally Toeplitz (GLT) sequences, respectively. In the context of structured matrices, the convergence properties of iterative methods, like multigrid or preconditioned Krylov techniques, are strictly related to the notion of symbol, a function whose role relies in describing the asymptotical distribution of the spectrum. This thesis can be seen as a byproduct of the combined use of powerful tools like symbol, spectral distribution, and GLT, when dealing with the numerical solution of structured linear systems. We approach such an issue both from a theoretical and practical viewpoint. On the one hand, we enlarge some known spectral distribution tools by proving the eigenvalue distribution of matrix-sequences obtained as combination of some algebraic operations on multilevel block Toeplitz matrices. On the other hand, we take advantage of the obtained results for designing efficient preconditioning techniques. Moreover, we focus on the numerical solution of structured linear systems coming from the following applications: image deblurring, fractional diffusion equations, and coupled PDEs. A spectral analysis of the arising structured sequences allows us either to study the convergence and predict the behavior of preconditioned Krylov and multigrid methods applied to the coefficient matrices, or to design effective preconditioners and multigrid solvers for the associated linear systems