Arvind K. Saibaba

Tania Bakhos, Peter Kitanidis, Scott Ladenheim et al.

An implementation of GMRES with multiple preconditioners (MPGMRES) is proposed for solving shifted linear systems with shift-and-invert preconditioners. With this type of preconditioner, the Krylov subspace can be built without requiring the matrix-vector product with the shifted matrix. Furthermore, the multipreconditioned search space is shown to grow only linearly with the number of preconditioners. This allows for a more efficient implementation of the algorithm. The proposed implementation is tested on shifted systems that arise in computational hydrology and the evaluation of different matrix functions. The numerical results indicate the effectiveness of the proposed approach.

Arvind K. Saibaba, Alen Alexanderian, Ilse C. F. Ipsen

We present randomized algorithms for estimating the trace and deter- minant of Hermitian positive semi-definite matrices. The algorithms are based on subspace iteration, and access the matrix only through matrix vector products. We analyse the error due to randomization, for starting guesses whose elements are Gaussian or Rademacher random variables. The analysis is cleanly separated into a structural (deterministic) part followed by a probabilistic part. Our absolute bounds for the expectation and concentration of the estimators are non-asymptotic and informative even for matrices of low dimension. For the trace estimators, we also present asymptotic bounds on the number of samples (columns of the starting guess) required to achieve a user-specified relative error. Numerical experiments illustrate the performance of the estimators and the tightness of the bounds on low-dimensional matrices; and on a challenging application in uncertainty quantification arising from Bayesian optimal experimental design.

Alen Alexanderian, Arvind K. Saibaba

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.

Zlatko Drmač, Arvind K. Saibaba

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.

Arvind K. Saibaba

This paper is concerned with the analysis of the randomized subspace iteration for the computation of low-rank approximations. We present three different kinds of bounds. First, we derive both bounds for the canonical angles between the exact and the approximate singular subspaces. Second, we derive bounds for the low-rank approximation in any unitarily invariant norm (including the Schatten-p norm). This generalizes the bounds for Spectral and Frobenius norms found in the literature. Third, we present bounds for the accuracy of the singular values. The bounds are structural in that they are applicable to any starting guess, be it random or deterministic, that satisfies some minimal assumptions. Specialized bounds are provided when a Gaussian random matrix is used as the starting guess. Numerical experiments demonstrate the effectiveness of the proposed bounds.

Arvind K. Saibaba, Johnathan Bardsley, D. Andrew Brown et al.

Hierarchical models in Bayesian inverse problems are characterized by an assumed prior probability distribution for the unknown state and measurement error precision, and hyper-priors for the prior parameters. Combining these probability models using Bayes' law often yields a posterior distribution that cannot be sampled from directly, even for a linear model with Gaussian measurement error and Gaussian prior. Gibbs sampling can be used to sample from the posterior, but problems arise when the dimension of the state is large. This is because the Gaussian sample required for each iteration can be prohibitively expensive to compute, and because the statistical efficiency of the Markov chain degrades as the dimension of the state increases. The latter problem can be mitigated using marginalization-based techniques, but these can be computationally prohibitive as well. In this paper, we combine the low-rank techniques of Brown, Saibaba, and Vallelian (2018) with the marginalization approach of Rue and Held (2005). We consider two variants of this approach: delayed acceptance and pseudo-marginalization. We provide a detailed analysis of the acceptance rates and computational costs associated with our proposed algorithms, and compare their performances on two numerical test cases---image deblurring and inverse heat equation.

Ege Ozsar, Misha Kilmer, Eric Miller et al.

We introduce PaLEnTIR, a significantly enhanced parametric level-set (PaLS) method addressing the restoration and reconstruction of piecewise constant objects. Our key contribution involves a unique PaLS formulation utilizing a single level-set function to restore scenes containing multi-contrast piecewise-constant objects without requiring knowledge of the number of objects or their contrasts. Unlike standard PaLS methods employing radial basis functions (RBFs), our model integrates anisotropic basis functions (ABFs), thereby expanding its capacity to represent a wider class of shapes. Furthermore, PaLEnTIR improves the conditioning of the Jacobian matrix, required as part of the parameter identification process, and consequently accelerates optimization methods. We validate PaLEnTIR's efficacy through diverse experiments encompassing sparse and limited angle of view X-ray computed tomography (2D and 3D), nonlinear diffuse optical tomography (DOT), denoising, and deconvolution tasks using both real and simulated data sets.

Ilse C. F. Ipsen, Arvind K. Saibaba

We consider a variety of criteria for selecting k representative columns from a real mxn matrix A, when sufficiently few columns are required, i.e., 1<= k<= min{rank(A), m/3}. The criteria include the following optimization problems: absolute volume and S-optimality maximization; norm, pseudo-inverse norm, and condition minimization number in the two-norm, Frobenius norm and Schatten p-norms for p>2; stable rank maximization; and the new criterion of relative volume maximization, which is inversely proportional to a power of the condition number. We show that these criteria are NP hard and many do not admit polynomial time approximation schemes (PTAS). To formulate the optimization problems as decision problems, we derive optimal values for the subset selection criteria, as well as expressions for partitioned pseudo-inverses. The results for minimization of the pseudo-inverse in the Frobenius norm are applicable to trace optimization in A-optimal design.

D. Andrew Brown, Arvind Saibaba, Sarah Vallélian

In Bayesian inverse problems, the posterior distribution is used to quantify uncertainty about the reconstructed solution. In practice, Markov chain Monte Carlo algorithms often are used to draw samples from the posterior distribution. However, implementations of such algorithms can be computationally expensive. We present a computationally efficient scheme for sampling high-dimensional Gaussian distributions in ill-posed Bayesian linear inverse problems. Our approach uses Metropolis-Hastings independence sampling with a proposal distribution based on a low-rank approximation of the prior-preconditioned Hessian. We show the dependence of the acceptance rate on the number of eigenvalues retained and discuss conditions under which the acceptance rate is high. We demonstrate our proposed sampler by using it with Metropolis-Hastings-within-Gibbs sampling in numerical experiments in image deblurring, computerized tomography, and NMR relaxometry.

Arvind K. Saibaba, Julianne Chung, Katrina Petroske

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.

Arvind K. Saibaba, Bhisham Dev Verma, Grey Ballard

Randomization has emerged as a powerful set of tools for large-scale matrix and tensor decompositions. Randomized algorithms involve computing sketches with random matrices. A prevalent approach is to take the random matrix as a standard Gaussian random matrix, for which the theory is well developed. However, this approach has the drawback that the cost of generating and multiplying by the random matrix can be prohibitively expensive. Khatri-Rao random projections (KRPs), obtained by sketching with Khatri-Rao products of random matrices, offer a viable alternative and are much cheaper to generate. However, the theoretical guarantees of using KRPs are much more pessimistic compared to their accuracy observed in practice. We attempt to close this gap by obtaining improved analysis of the use of KRPs in matrix and tensor low-rank decompositions. We propose and analyze a new algorithm for low-rank approximations of block-structured matrices (e.g., block Hankel) using KRPs. We also show how to accelerate tensor computations in the Tucker format using KRPs and give theoretical guarantees of the resulting low-rank approximations. Numerical experiments on synthetic and real-world tensors show the computational benefits of the proposed methods.

Arvind K. Saibaba, Ilse C. F. Ipsen

We consider matrices $\boldsymbol{A}(\boldsymbolθ)\in\mathbb{R}^{m\times m}$ that depend, possibly nonlinearly, on a parameter $\boldsymbolθ$ from a compact parameter space $Θ$. We present a Monte Carlo estimator for minimizing $\text{trace}(\boldsymbol{A}(\boldsymbolθ))$ over all $\boldsymbolθ\inΘ$, and determine the sampling amount so that the backward error of the estimator is bounded with high probability. We derive two types of bounds, based on epsilon nets and on generic chaining. Both types predict a small sampling amount for matrices $\boldsymbol{A}(\boldsymbolθ)$ with small offdiagonal mass, and parameter spaces $Θ$ of small ``size.'' Dependence on the matrix dimension~$m$ is only weak or not explicit. The bounds based on epsilon nets are easier to evaluate and come with fully specified constants. In contrast, the bounds based on chaining depend on the Talagrand functionals which are difficult to evaluate, except in very special cases. Comparisons between the two types of bounds are difficult, although the literature suggests that chaining bounds can be superior.

Eric Chi, Liuiyi Hu, Arvind K. Saibaba et al.

We consider the problem of performing matrix completion with side information on row-by-row and column-by-column similarities. We build upon recent proposals for matrix estimation with smoothness constraints with respect to row and column graphs. We present a novel iterative procedure for directly minimizing an information criterion in order to select an appropriate amount row and column smoothing, namely perform model selection. We also discuss how to exploit the special structure of the problem to scale up the estimation and model selection procedure via the Hutchinson estimator. We present simulation results and an application to predicting associations in imaging-genomics studies.

Arvind K. Saibaba, Jonghyun Lee, Peter K. Kitanidis

We describe randomized algorithms for computing the dominant eigenmodes of the Generalized Hermitian Eigenvalue Problem (GHEP) $Ax=λBx$, with $A$ Hermitian and $B$ Hermitian and positive definite. The algorithms we describe only require forming operations $Ax$, $Bx$ and $B^{-1}x$ and avoid forming square-roots of $B$ (or operations of the form, $B^{1/2}x$ or $B^{-1/2}x$). We provide a convergence analysis and a posteriori error bounds that build upon the work of~\cite{halko2011finding,liberty2007randomized,martinsson2011randomized} (which have been derived for the case $B=I$). Additionally, we derive some new results that provide insight into the accuracy of the eigenvalue calculations. The error analysis shows that the randomized algorithm is most accurate when the generalized singular values of $B^{-1}A$ decay rapidly. A randomized algorithm for the Generalized Singular Value Decomposition (GSVD) is also provided. Finally, we demonstrate the performance of our algorithm on computing the Karhunen-Loève expansion, which is a computationally intensive GHEP problem with rapidly decaying eigenvalues.

Arvind K. Saibaba, Misha Kilmer, Eric Miller et al.

The image reconstruction of chromophore concentrations using Diffuse Optical Tomography (DOT) data can be described mathematically as an ill-posed inverse problem. Recent work has shown that the use of hyperspectral DOT data, as opposed to data sets comprising of a single or, at most, a dozen wavelengths, has the potential for improving the quality of the reconstructions. The use of hyperspectral diffuse optical data in the formulation and solution of the inverse problem poses a significant computational burden. The forward operator is, in actuality, nonlinear. However, under certain assumptions, a linear approximation, called the Born approximation, provides a suitable surrogate for the forward operator, and we assume this to be true in the present work. Computation of the Born matrix requires the solution of thousands of large scale discrete PDEs and the reconstruction problem, requires matrix-vector products with the (dense) Born matrix. In this paper, we address both of these difficulties, thus making the Born approach a computational viable approach for hyDOT reconstruction. In this paper, we assume that the images we wish to reconstruct are anomalies of unknown shape and constant value, described using a parametric level set approach, (PaLS) on a constant background. Specifically, to address the issue of the PDE solves, we develop a novel recycling-based Krylov subspace approach that leverages certain system similarities across wavelengths. To address expense of using the Born operator in the inversion, we present a fast algorithm for compressing the Born operator that locally compresses across wavelengths for a given source-detector set and then recursively combines the low-rank factors to provide a global low-rank approximation. This low-rank approximation can be used implicitly to speed up the recovery of the shape parameters and the chromophore concentrations.