A Householder-based algorithm for Hessenberg-triangular reduction

The QZ algorithm for computing eigenvalues and eigenvectors of a matrix pencil $A - \lambda B$ requires that the matrices first be reduced to Hessenberg-triangular (HT) form. The current method of choice for HT reduction relies entirely on Givens rotations regrouped and accumulated into small dense matrices which are subsequently applied using matrix multiplication routines. A non-vanishing fraction of the total flop count must nevertheless still be performed as sequences of overlapping Givens rotations alternately applied from the left and from the right. The many data dependencies associated with this computational pattern leads to inefficient use of the processor and poor scalability. In this paper, we therefore introduce a fundamentally different approach that relies entirely on (large) Householder reflectors partially accumulated into block reflectors, by using (compact) WY representations. Even though the new algorithm requires more floating point operations than the state of the art algorithm, extensive experiments on both real and synthetic data indicate that it is still competitive, even in a sequential setting. The new algorithm is conjectured to have better parallel scalability, an idea which is partially supported by early small-scale experiments using multi-threaded BLAS. The design and evaluation of a parallel formulation is future work

MATHICSE Technical Report : A continuation multilevel Monte Carlo algorithm

We propose a novel Continuation Multi Level Monte Carlo (CMLMC) algorithm for weak approximation of stochastic models that are described in terms of differential equations either driven by random measures or with random coefficients. The CMLMC algorithm solves the given approximation problem for a sequence of decreasing tolerances, ending with the desired one. CMLMC assumes discretization hierarchies that are defined a priori for each level and are geometrically refined across levels. The actual choice of computational work across levels is based on parametric models for the average cost per sample and the corresponding weak and strong errors. These parameters are calibrated using Bayesian estimation, taking particular notice of the deepest levels of the discretization hierarchy, where only few realizations are available to produce the estimates. The resulting CMLMC estimator exhibits a nontrivial splitting between bias and statistical contributions. We also show the asymptotic normality of the statistical error in the MLMC estimator and justify in this way our error estimate that allows prescribing both required accuracy and confidence in the final result. Numerical examples substantiate the above results and illustrate the corresponding computational savings

MATHICSE Technical Report : Optimization of mesh hierarchies in multilevel Monte Carlo samplers

We perform a general optimization of the parameters in the Multilevel Monte Carlo (MLMC) discretization hierarchy based on uniform discretization methods with general approximation orders and computational costs. Moreover, we discuss extensions to non-uniform discretizations based on a priori renements and the effect of imposing constraints on the largest and/or smallest mesh sizes. We optimize geometric and non-geometric hierarchies and compare them to each other, concluding that the geometric hierarchies, when optimized, are nearly optimal and have the same asymptotic computational complexity. We discuss how enforcing domain constraints on parameters of MLMC hierarchies affects the opti- mality of these hierarchies. These domain constraints include an upper and lower bound on the mesh size or enforcing that the number of samples and the number of discretization elements are integers. We also discuss the optimal tolerance splitting between the bias and the statistical error contributions and its asymp- totic behavior. To provide numerical grounds for our theoretical results, we apply these optimized hierarchies together with the Continuation MLMC Algorithm [13] that we recently developed, to several examples. These include the approxima- tion of three-dimensional elliptic partial differential equations with random inputs based on FEM with either direct or iterative solvers and It^o stochastic differential equations based on the Milstein scheme

MATHICSE Technical Report : Convergence of quasi-optimal sparse grid approximation of Hilbert-valued functions: application to random elliptic PDEs

In this work we provide a convergence analysis for the quasi-optimal version of the Stochastic Sparse Grid Collocation method we had presented in our previous work \On the optimal polynomial approximation of Stochastic PDEs by Galerkin and Collocation methods" [6]. Here the construction of a sparse grid is recast into a knapsack problem: a profit is assigned to each hi- erarchical surplus and only the most profitable ones are added to the sparse grid. The convergence rate of the sparse grid approximation error with respect to the number of points in the grid is then shown to depend on weighted summability properties of the sequence of profits. This argument is very gen- eral and can be applied to sparse grids built with any uni-variate family of points, both nested and non-nested. As an example, we apply such quasi-optimal sparse grid to the solution of a particular elliptic PDE with stochastic diffusion coefficients, namely the \inclusions problem": we detail the conver- gence estimate obtained in this case, using polynomial interpolation on either nested (Clenshaw{Curtis) or non-nested (Gauss{Legendre) abscissas, verify its sharpness numerically, and compare the performance of the resulting quasi- optimal grids with a few alternative sparse grids construction schemes recently proposed in literature

A Parallel QZ Algorithm For Distributed Memory HPC Systems

Appearing frequently in applications, generalized eigenvalue problems represent one of the core problems in numerical linear algebra. The QZ algorithm of Moler and Stewart is the most widely used algorithm for addressing such problems. Despite its importance, little attention has been paid to the parallelization of the QZ algorithm. The purpose of this work is to fill this gap. We propose a parallelization of the QZ algorithm that incorporates all modern ingredients of dense eigensolvers, such as multishift and aggressive early deflation techniques. To deal with (possibly many) infinite eigenvalues, a new parallel deflation strategy is developed. Numerical experiments for several random and application examples demonstrate the effectiveness of our algorithm on two different distributed memory HPC systems