A comparative study of null-space factorizations for sparse symmetric saddle point systems

Null-space methods for solving saddle point systems of equations have long been used to transform an indefinite system into a symmetric positive definite one of smaller dimension. A number of independent works in the literature have identified that we can interpret a null-space method as a matrix factorization. We review these findings, highlight links between them, and bring them into a unified framework. We also investigate the suitability of using null-space factorizations to derive sparse direct methods, and present numerical results for both practical and academic problems

Implicit QR algorithms for palindromic and even eigenvalue problems

In the spirit of the Hamiltonian QR algorithm and other bidirectional chasing algorithms, a structure-preserving variant of the implicit QR algorithm for palindromic eigenvalue problems is proposed. This new palindromic QR algorithm is strongly backward stable and requires less operations than the standard QZ algorithm, but is restricted to matrix classes where a preliminary reduction to structured Hessenberg form can be performed. By an extension of the implicit Q theorem, the palindromic QR algorithm is shown to be equivalent to a previously developed explicit version. Also, the classical convergence theory for the QR algorithm can be extended to prove local quadratic convergence. We briefly demonstrate how even eigenvalue problems can be addressed by similar techniques. © 2008 Springer Science+Business Media, LLC

A Kogbetliantz-type algorithm for the hyperbolic SVD

In this paper a two-sided, parallel Kogbetliantz-type algorithm for the hyperbolic singular value decomposition (HSVD) of real and complex square matrices is developed, with a single assumption that the input matrix, of order $n$, admits such a decomposition into the product of a unitary, a non-negative diagonal, and a $J$-unitary matrix, where $J$ is a given diagonal matrix of positive and negative signs. When $J=\pm I$, the proposed algorithm computes the ordinary SVD. The paper's most important contribution -- a derivation of formulas for the HSVD of $2\times 2$ matrices -- is presented first, followed by the details of their implementation in floating-point arithmetic. Next, the effects of the hyperbolic transformations on the columns of the iteration matrix are discussed. These effects then guide a redesign of the dynamic pivot ordering, being already a well-established pivot strategy for the ordinary Kogbetliantz algorithm, for the general, $n\times n$ HSVD. A heuristic but sound convergence criterion is then proposed, which contributes to high accuracy demonstrated in the numerical testing results. Such a $J$-Kogbetliantz algorithm as presented here is intrinsically slow, but is nevertheless usable for matrices of small orders.Comment: a heavily revised version with 32 pages and 4 figure

MATHICSE Technical Report : Convergence estimates in probability and in expectation for discrete least squares with noisy evaluations at random points

We study the accuracy of the discrete least-squares approximation on a finite dimensional space of a real-valued target function from noisy pointwise evaluations at independent random points distributed according to a given sampling probability measure. The convergence estimates are given in mean-square sense with respect to the sampling measure. The noise may be correlated with the location of the evaluation and may have nonzero mean (offset). We consider both cases of bounded or square-integrable noise / offset. We prove conditions between the number of sampling points and the dimension of the underlying approximation space that ensure a stable and accurate approximation. Particular focus is on deriving estimates in probability within a given confidence level. We analyze how the best approximation error and the noise terms affect the convergence rate and the overall confidence level achieved by the convergence estimate. The proofs of our convergence estimates in probability use arguments from the theory of large deviations to bound the noise term. Finally we address the particular case of multivariate polynomial approximation spaces with any density in the beta family, including uniform and Chebyshev

MATHICSE Technical Report : An adaptive sparse grid algorithm for elliptic PDEs with lognormal diffusion coefficient

In this work we build on the classical adaptive sparse grid algorithm (T. Gerstner and M. Griebel, Dimension-adaptive tensor-product quadrature), obtaining an enhanced version capable of using non-nested collocation points, and supporting quadrature and interpolation on unbounded sets. We also consider several profit indicators that are suitable to drive the adaptation process.We then use such algorithm to solve an important test case in Uncertainty Quantification problem, namely the Darcy equation with lognormal permeability random field, and compare the results with those obtained with the quasi-optimal sparse grids based on profit estimates, which we have proposed in our previous works (cf. e.g. Convergence of quasi-optimal sparse grids approximation of Hilbert-valued functions: application to random elliptic PDEs). To treat the case of rough permeability fields, in which a sparse grid approach may not be suitable, we propose to use the adaptive sparse grid quadrature as a control variate in a Monte Carlo simulation. Numerical results show that the adaptive sparse grids have performances similar to those of the quasi-optimal sparse grids and are very effective in the case of smooth permeability fields. Moreover, their use as control variate in a Monte Carlo simulation allows to tackle efficiently also problems with rough coefficients, significantly improving the performances of a standard Monte Carlo scheme

Accurate solution of Bayesian inverse uncertainty quantification problems combining reduced basis methods and reduction error models

Computational inverse problems related to partial differential equations (PDEs) often contain nuisance parameters that cannot be effectively identified but still need to be considered as part of the problem. The objective of this work is to show how to take advantage of a reduced order framework to speed up Bayesian inversion on the identifiable parameters of the system, while marginalizing away the (potentially large number of) nuisance parameters. The key ingredients are twofold. On the one hand, we rely on a reduced basis (RB) method, equipped with computable a posteriori error bounds, to speed up the solution of the forward problem. On the other hand, we develop suitable reduction error models (REMs) to quantify in an inexpensive way the error between the full-order and the reduced-order approximation of the forward problem, in order to gauge the effect of this error on the posterior distribution of the identifiable parameters. Numerical results dealing with inverse problems governed by elliptic PDEs in the case of both scalar parameters and parametric fields highlight the combined role played by RB accuracy and REM effectivity

MATHICSE Technical Report : A multi level Monte Carlo method with control variate for elliptic PDEs with log-normal coefficients

We consider the numerical approximation of the stochastic Darcy problem with log-normal permeability field and propose a novel Multi Level Monte Carlo approach with a control variate variance reduction technique on each level. We model the log-permeability as a stationary Gaussian random field with a covariance function belonging to the so called Matérn family, which includes both fields with very limited and very high spatial regularity. The control variate is obtained starting from the solution of an auxiliary problem with smoothed permeability coefficient and its expected value is effectively computed with a Stochastic Collocation method on the finest level in which the control variate is applied. We analyze the variance reduction induced by the control variate, and the total mean square error of the new estimator. To conclude we present some numerical examples and a comparison with the standard Multi Level Monte Carlo method, which shows the effectiveness of the proposed method

A block algorithm for computing antitriangular factorizations of symmetric matrices

Any symmetric matrix can be reduced to antitriangular form in finitely many steps by orthogonal similarity transformations. This form reveals the inertia of the matrix and has found applications in, e.g., model predictive control and constraint preconditioning. Originally proposed by Mastronardi and Van Dooren, the existing algorithm for performing the reduction to antitriangular form is primarily based on Householder reflectors and Givens rotations. The poor memory access pattern of these operations implies that the performance of the algorithm is bound by the memory bandwidth. In this work, we develop a block algorithm that performs all operations almost entirely in terms of level 3 BLAS operations, which feature a more favorable memory access pattern and lead to better performance. These performance gains are confirmed by numerical experiments that cover a wide range of different inertia