HOID: Higher Order Interpolatory Decomposition for tensors based on Tucker representation

We derive a CUR-type factorization for tensors in the Tucker format based on interpolatory decomposition, which we will denote as Higher Order Interpolatory Decomposition (HOID). Given a tensor $\mathcal{X}$, the algorithm provides a set of column vectors $\{ \mathbf{C}_n\}_{n=1}^d$ which are columns extracted from the mode-$n$ tensor unfolding, along with a core tensor $\mathcal{G}$ and together, they satisfy some error bounds. Compared to the Higher Order SVD (HOSVD) algorithm, the HOID provides a decomposition that preserves certain important features of the original tensor such as sparsity, non-negativity, integer values, etc. Error bounds along with detailed estimates of computational costs are provided. The algorithms proposed in this paper have been validated against carefully chosen numerical examples which highlight the favorable properties of the algorithms. Related methods for subset selection proposed for matrix CUR decomposition, such as Discrete Empirical Interpolation method (DEIM) and leverage score sampling, have also been extended to tensors and are compared against our proposed algorithms.Comment: 28 pages, 9 figures, minor revision

Generalized hybrid iterative methods for large-scale Bayesian inverse problems

We develop a generalized hybrid iterative approach for computing solutions to large-scale Bayesian inverse problems. We consider a hybrid algorithm based on the generalized Golub-Kahan bidiagonalization for computing Tikhonov regularized solutions to problems where explicit computation of the square root and inverse of the covariance kernel for the prior covariance matrix is not feasible. This is useful for large-scale problems where covariance kernels are defined on irregular grids or are only available via matrix-vector multiplication, e.g., those from the Mat\'{e}rn class. We show that iterates are equivalent to LSQR iterates applied to a directly regularized Tikhonov problem, after a transformation of variables, and we provide connections to a generalized singular value decomposition filtered solution. Our approach shares many benefits of standard hybrid methods such as avoiding semi-convergence and automatically estimating the regularization parameter. Numerical examples from image processing demonstrate the effectiveness of the described approaches

Efficient D-optimal design of experiments for infinite-dimensional Bayesian linear inverse problems

We develop a computational framework for D-optimal experimental design for PDE-based Bayesian linear inverse problems with infinite-dimensional parameters. We follow a formulation of the experimental design problem that remains valid in the infinite-dimensional limit. The optimal design is obtained by solving an optimization problem that involves repeated evaluation of the log-determinant of high-dimensional operators along with their derivatives. Forming and manipulating these operators is computationally prohibitive for large-scale problems. Our methods exploit the low-rank structure in the inverse problem in three different ways, yielding efficient algorithms. Our main approach is to use randomized estimators for computing the D-optimal criterion, its derivative, as well as the Kullback--Leibler divergence from posterior to prior. Two other alternatives are proposed based on a low-rank approximation of the prior-preconditioned data misfit Hessian, and a fixed low-rank approximation of the prior-preconditioned forward operator. Detailed error analysis is provided for each of the methods, and their effectiveness is demonstrated on a model sensor placement problem for initial state reconstruction in a time-dependent advection-diffusion equation in two space dimensions.Comment: 27 pages, 9 figure

The Discrete Empirical Interpolation Method: Canonical Structure and Formulation in Weighted Inner Product Spaces

New contributions are offered to the theory and practice of the Discrete Empirical Interpolation Method (DEIM). These include a detailed characterization of the canonical structure; a substantial tightening of the error bound for the DEIM oblique projection, based on index selection via a strong rank revealing QR factorization; and an extension of the DEIM approximation to weighted inner products defined by a real symmetric positive-definite matrix $W$. The weighted DEIM ($W$-DEIM) can be deployed in the more general framework where the POD Galerkin projection is formulated in a discretization of a suitable energy inner product such that the Galerkin projection preserves important physical properties such as e.g. stability. Also, a special case of $W$-DEIM is introduced, which is DGEIM, a discrete version of the Generalized Empirical Interpolation Method that allows generalization of the interpolation via a dictionary of linear functionals.Comment: 26 pages, 6 figure

Randomized Discrete Empirical Interpolation Method for Nonlinear Model Reduction

Discrete empirical interpolation method (DEIM) is a popular technique for nonlinear model reduction and it has two main ingredients: an interpolating basis that is computed from a collection of snapshots of the solution and a set of indices which determine the nonlinear components to be simulated. The computation of these two ingredients dominates the overall cost of the DEIM algorithm. To specifically address these two issues, we present randomized versions of the DEIM algorithm. There are three main contributions of this paper. First, we use randomized range finding algorithms to efficiently find an approximate DEIM basis. Second, we develop randomized subset selection tools, based on leverage scores, to efficiently select the nonlinear components. Third, we develop several theoretical results that quantify the accuracy of the randomization on the DEIM approximation. We also present numerical experiments that demonstrate the benefits of the proposed algorithms.Comment: 27 pages, 8 figure

Uncertainty quantification in large Bayesian linear inverse problems using Krylov subspace methods

For linear inverse problems with a large number of unknown parameters, uncertainty quantification remains a challenging task. In this work, we use Krylov subspace methods to approximate the posterior covariance matrix and describe efficient methods for exploring the posterior distribution. Assuming that Krylov methods (e.g., based on the generalized Golub-Kahan bidiagonalization) have been used to compute an estimate of the solution, we get an approximation of the posterior covariance matrix for `free.' We provide theoretical results that quantify the accuracy of the approximation and of the resulting posterior distribution. Then, we describe efficient methods that use the approximation to compute measures of uncertainty, including the Kullback-Liebler divergence. We present two methods that use preconditioned Lanczos methods to efficiently generate samples from the posterior distribution. Numerical examples from tomography demonstrate the effectiveness of the described approaches.Comment: 26 pages, 4 figures, 2 tables. Under revie

Monte Carlo Estimators for the Schatten p-norm of Symmetric Positive Semidefinite Matrices

We present numerical methods for computing the Schatten $p$-norm of positive semi-definite matrices. Our motivation stems from uncertainty quantification and optimal experimental design for inverse problems, where the Schatten $p$-norm defines a design criterion known as the P-optimal criterion. Computing the Schatten $p$-norm of high-dimensional matrices is computationally expensive. We propose a matrix-free method to estimate the Schatten $p$-norm using a Monte Carlo estimator and derive convergence results and error estimates for the estimator. To efficiently compute the Schatten $p$-norm for non-integer and large values of $p$, we use an estimator using a Chebyshev polynomial approximation and extend our convergence and error analysis to this setting as well. We demonstrate the performance of our proposed estimators on several test matrices and through an application to optimal experimental design of a model inverse problem.Comment: 21 pages, 10 figures, 1 tabl

Goal-Oriented Optimal Design of Experiments for Large-Scale Bayesian Linear Inverse Problems

We develop a framework for goal-oriented optimal design of experiments (GOODE) for large-scale Bayesian linear inverse problems governed by PDEs. This framework differs from classical Bayesian optimal design of experiments (ODE) in the following sense: we seek experimental designs that minimize the posterior uncertainty in the experiment end-goal, e.g., a quantity of interest (QoI), rather than the estimated parameter itself. This is suitable for scenarios in which the solution of an inverse problem is an intermediate step and the estimated parameter is then used to compute a QoI. In such problems, a GOODE approach has two benefits: the designs can avoid wastage of experimental resources by a targeted collection of data, and the resulting design criteria are computationally easier to evaluate due to the often low-dimensionality of the QoIs. We present two modified design criteria, A-GOODE and D-GOODE, which are natural analogues of classical Bayesian A- and D-optimal criteria. We analyze the connections to other ODE criteria, and provide interpretations for the GOODE criteria by using tools from information theory. Then, we develop an efficient gradient-based optimization framework for solving the GOODE optimization problems. Additionally, we present comprehensive numerical experiments testing the various aspects of the presented approach. The driving application is the optimal placement of sensors to identify the source of contaminants in a diffusion and transport problem. We enforce sparsity of the sensor placements using an $\ell_1$-norm penalty approach, and propose a practical strategy for specifying the associated penalty parameter.Comment: 25 pages, 13 figure

Approximating monomials using Chebyshev polynomials

This paper considers the approximation of a monomial $x^n$ over the interval $[-1,1]$ by a lower-degree polynomial. This polynomial approximation can be easily computed analytically and is obtained by truncating the analytical Chebyshev series expansion of $x^n$. The error in the polynomial approximation in the supremum norm has an exact expression with an interesting probabilistic interpretation. We use this interpretation along with concentration inequalities to develop a useful upper bound for the error.Comment: 6 pages, 2 figure

Randomized algorithms for low-rank tensor decompositions in the Tucker format

Many applications in data science and scientific computing involve large-scale datasets that are expensive to store and compute with, but can be efficiently compressed and stored in an appropriate tensor format. In recent years, randomized matrix methods have been used to efficiently and accurately compute low-rank matrix decompositions. Motivated by this success, we focus on developing randomized algorithms for tensor decompositions in the Tucker representation. Specifically, we present randomized versions of two well-known compression algorithms, namely, HOSVD and STHOSVD. We present a detailed probabilistic analysis of the error of the randomized tensor algorithms. We also develop variants of these algorithms that tackle specific challenges posed by large-scale datasets. The first variant adaptively finds a low-rank representation satisfying a given tolerance and it is beneficial when the target-rank is not known in advance. The second variant preserves the structure of the original tensor, and is beneficial for large sparse tensors that are difficult to load in memory. We consider several different datasets for our numerical experiments: synthetic test tensors and realistic applications such as the compression of facial image samples in the Olivetti database and word counts in the Enron email dataset.Comment: 28 pages, 4 figures, 6 table