Source author record

Arvind K. Saibaba

Arvind K. Saibaba appears in the imported research catalog. Authorship, coauthor and topic links are available while profile ownership is still unclaimed.

ResearcherUnclaimed source record

Catalog footprint

What is connected

16works
3topics
4close collaborators

Actions

Connect this record

Log in to claim

Research graph

See the researcher in context

Open full explorer

Inspect adjacent papers, topics, institutions and collaborators without losing the researcher page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

Published work

16 published item(s)

preprint2022arXiv

Efficient algorithms for Bayesian Inverse Problems with Whittle--Matérn Priors

This paper tackles efficient methods for Bayesian inverse problems with priors based on Whittle--Matérn Gaussian random fields. The Whittle--Matérn prior is characterized by a mean function and a covariance operator that is taken as a negative power of an elliptic differential operator. This approach is flexible in that it can incorporate a wide range of prior information including non-stationary effects, but it is currently computationally advantageous only for integer values of the exponent. In this paper, we derive an efficient method for handling all admissible noninteger values of the exponent. The method first discretizes the covariance operator using finite elements and quadrature, and uses preconditioned Krylov subspace solvers for shifted linear systems to efficiently apply the resulting covariance matrix to a vector. This approach can be used for generating samples from the distribution in two different ways: by solving a stochastic partial differential equation, and by using a truncated Karhunen-Loève expansion. We show how to incorporate this prior representation into the infinite-dimensional Bayesian formulation, and show how to efficiently compute the maximum a posteriori estimate, and approximate the posterior variance. Although the focus of this paper is on Bayesian inverse problems, the techniques developed here are applicable to solving systems with fractional Laplacians and Gaussian random fields. Numerical experiments demonstrate the performance and scalability of the solvers and their applicability to model and real-data inverse problems in tomography and a time-dependent heat equation.

preprint2022arXiv

Robust Parameter Identifiability Analysis via Column Subset Selection

We advocate a numerically reliable and accurate approach for practical parameter identifiability analysis: Applying column subset selection (CSS) to the sensitivity matrix, instead of computing an eigenvalue decomposition of the Fischer information matrix. Identifiability analysis via CSS has three advantages: (i) It quantifies reliability of the subsets of parameters selected as identifiable and unidentifiable. (ii) It establishes criteria for comparing the accuracy of different algorithms. (iii) The implementations are numerically more accurate and reliable than eigenvalue methods applied to the Fischer matrix, yet without an increase in computational cost. The effectiveness of the CSS methods is illustrated with extensive numerical experiments on sensitivity matrices from six physical models, as well as on adversarial synthetic matrices. Among the CSS methods, we recommend an implementation based on the strong rank-revealing QR algorithm because of its rigorous accuracy guarantees for both identifiable and non-identifiable parameters.

preprint2021arXiv

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.

preprint2020arXiv

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.

preprint2020arXiv

Randomization and reweighted $\ell_1$-minimization for A-optimal design of linear inverse problems

We consider optimal design of PDE-based Bayesian linear inverse problems with infinite-dimensional parameters. We focus on the A-optimal design criterion, defined as the average posterior variance and quantified by the trace of the posterior covariance operator. We propose using structure exploiting randomized methods to compute the A-optimal objective function and its gradient, and provide a detailed analysis of the error for the proposed estimators. To ensure sparse and binary design vectors, we develop a novel reweighted $\ell_1$-minimization algorithm. We also introduce a modified A-optimal criterion and present randomized estimators for its efficient computation. We present numerical results illustrating the proposed methods on a model contaminant source identification problem, where the inverse problem seeks to recover the initial state of a contaminant plume, using discrete measurements of the contaminant in space and time.

preprint2020arXiv

Randomized Algorithms for Generalized Singular Value Decomposition with Application to Sensitivity Analysis

The generalized singular value decomposition (GSVD) is a valuable tool that has many applications in computational science. However, computing the GSVD for large-scale problems is challenging. Motivated by applications in hyper-differential sensitivity analysis (HDSA), we propose new randomized algorithms for computing the GSVD which use randomized subspace iteration and weighted QR factorization. Detailed error analysis is given which provides insight into the accuracy of the algorithms and the choice of the algorithmic parameters. We demonstrate the performance of our algorithms on test matrices and a large-scale model problem where HDSA is used to study subsurface flow.

preprint2020arXiv

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.

preprint2016arXiv

A Randomized Tensor Singular Value Decomposition based on the t-product

The tensor Singular Value Decomposition (t-SVD) for third order tensors that was proposed by Kilmer and Martin~\cite{2011kilmer} has been applied successfully in many fields, such as computed tomography, facial recognition, and video completion. In this paper, we propose a method that extends a well-known randomized matrix method to the t-SVD. This method can produce a factorization with similar properties to the t-SVD, but is more computationally efficient on very large datasets. We present details of the algorithm, theoretical results, and provide numerical results that show the promise of our approach for compressing and analyzing datasets. We also present an improved analysis of the randomized subspace iteration for matrices, which may be of independent interest to the scientific community.

preprint2016arXiv

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.

preprint2016arXiv

Multipreconditioned GMRES for Shifted Systems

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.

preprint2015arXiv

A Fast Algorithm for Parabolic PDE-based Inverse Problems Based on Laplace Transforms and Flexible Krylov Solvers

We consider the problem of estimating parameters in large-scale weakly nonlinear inverse problems for which the underlying governing equations is a linear, time-dependent, parabolic partial differential equation. A major challenge in solving these inverse problems using Newton-type methods is the computational cost associated with solving the forward problem and with repeated construction of the Jacobian, which represents the sensitivity of the measurements to the unknown parameters. Forming the Jacobian can be prohibitively expensive because it requires repeated solutions of the forward and adjoint time-dependent parabolic partial differential equations corresponding to multiple sources and receivers. We propose an efficient method based on a Laplace transform-based exponential time integrator combined with a flexible Krylov subspace approach to solve the resulting shifted systems of equations efficiently. Our proposed solver speeds up the computation of the forward and adjoint problems, thus yielding significant speedup in total inversion time. We consider an application from Transient Hydraulic Tomography (THT), which is an imaging technique to estimate hydraulic parameters related to the subsurface from pressure measurements obtained by a series of pumping tests. The algorithms discussed are applied to a synthetic example taken from THT to demonstrate the resulting computational gains of this proposed method.

preprint2015arXiv

Fast computation of uncertainty quantification measures in the geostatistical approach to solve inverse problems

We consider the computational challenges associated with uncertainty quantification involved in parameter estimation such as seismic slowness and hydraulic transmissivity fields. The reconstruction of these parameters can be mathematically described as Inverse Problems which we tackle using the Geostatistical approach. The quantification of uncertainty in the Geostatistical approach involves computing the posterior covariance matrix which is prohibitively expensive to fully compute and store. We consider an efficient representation of the posterior covariance matrix at the maximum a posteriori (MAP) point as the sum of the prior covariance matrix and a low-rank update that contains information from the dominant generalized eigenmodes of the data misfit part of the Hessian and the inverse covariance matrix. The rank of the low-rank update is typically independent of the dimension of the unknown parameter. The cost of our method scales as $\bigO(m\log m)$ where $m $ dimension of unknown parameter vector space. Furthermore, we show how to efficiently compute measures of uncertainty that are based on scalar functions of the posterior covariance matrix. The performance of our algorithms is demonstrated by application to model problems in synthetic travel-time tomography and steady-state hydraulic tomography. We explore the accuracy of the posterior covariance on different experimental parameters and show that the cost of approximating the posterior covariance matrix depends on the problem size and is not sensitive to other experimental parameters.

preprint2015arXiv

Fast Kalman Filter using Hierarchical-matrices and low-rank perturbative approach

We develop a fast algorithm for Kalman Filter applied to the random walk forecast model. The key idea is an efficient representation of the estimate covariance matrix at each time-step as a weighted sum of two contributions - the process noise covariance matrix and a low rank term computed from a generalized eigenvalue problem, which combines information from the noise covariance matrix and the data. We describe an efficient algorithm to update the weights of the above terms and the computation of eigenmodes of the generalized eigenvalue problem (GEP). The resulting algorithm for the Kalman filter with a random walk forecast model scales as $\bigO(N)$ in memory and $\bigO(N \log N)$ in computational cost, where $N$ is the number of grid points. We show how to efficiently compute measures of uncertainty and conditional realizations from the state distribution at each time step. An extension to the case with nonlinear measurement operators is also discussed. Numerical experiments demonstrate the performance of our algorithms, which are applied to a synthetic example from monitoring CO$_2$ in the subsurface using travel time tomography.

preprint2015arXiv

Randomized algorithms for Generalized Hermitian Eigenvalue Problems with application to computing Karhunen-Loève expansion

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.

preprint2014arXiv

A Flexible Krylov Solver for Shifted Systems with Application to Oscillatory Hydraulic Tomography

We discuss efficient solutions to systems of shifted linear systems arising in computations for oscillatory hydraulic tomography (OHT). The reconstruction of hydrogeological parameters such as hydraulic conductivity and specific storage using limited discrete measurements of pressure (head) obtained from sequential oscillatory pumping tests, leads to a nonlinear inverse problem. We tackle this using the quasi-linear geostatistical approach \cite{kitanidis1995quasi}. This method requires repeated solution of the forward (and adjoint) problem for multiple frequencies, for which we use flexible preconditioned Krylov subspace solvers specifically designed for shifted systems based on ideas in \cite{gu2007flexible}. The solvers allow the preconditioner to change at each iteration. We analyze the convergence of the solver and perform an error analysis when an iterative solver is used for inverting the preconditioner matrices. Finally, we apply our algorithm to a challenging application taken from oscillatory hydraulic tomography to demonstrate the computational gains by using the resulting method.

preprint2014arXiv

Fast algorithms for hyperspectral Diffuse Optical Tomography

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.