@article{pasha_saibaba_gazzola_espanol_de sturler_2023, title={A COMPUTATIONAL FRAMEWORK FOR EDGE-PRESERVING REGULARIZATION IN DYNAMIC INVERSE PROBLEMS}, volume={58}, ISSN={["1068-9613"]}, DOI={10.1553/etna_vol58s486}, abstractNote={. We devise efficient methods for dynamic inverse problems, where both the quantities of interest and the forward operator (measurement process) may change in time. Our goal is to solve for all the quantities of interest simultaneously. We consider large-scale ill-posed problems made more challenging by their dynamic nature and, possibly, by the limited amount of available data per measurement step. To alleviate these difficulties, we apply a unified class of regularization methods that enforce simultaneous regularization in space and time (such as edge enhancement at each time instant and proximity at consecutive time instants) and achieve this with low computational cost and enhanced accuracy. More precisely, we develop iterative methods based on a majorization-minimization (MM) strategy with quadratic tangent majorant, which allows the resulting least-squares problem with a total variation regularization term to be solved with a generalized Krylov subspace (GKS) method; the regularization parameter can be determined automatically and efficiently at each iteration. Numerical examples from a wide range of applications, such as limited-angle computerized tomography (CT), space-time image deblurring, and photoacoustic tomography (PAT), illustrate the effectiveness of the described approaches.}, journal={ELECTRONIC TRANSACTIONS ON NUMERICAL ANALYSIS}, author={Pasha, Mirjeta and Saibaba, Arvind K. and Gazzola, Silvia and Espanol, Malena I. and De Sturler, Eric}, year={2023}, pages={486–516} }
 @article{hallman_ipsen_saibaba_2023, title={MONTE CARLO METHODS FOR ESTIMATING THE DIAGONAL OF A REAL SYMMETRIC MATRIX}, volume={44}, ISSN={["1095-7162"]}, url={https://doi.org/10.1137/22M1476277}, DOI={10.1137/22M1476277}, abstractNote={For real symmetric matrices that are accessible only through matrix vector products, we present Monte Carlo estimators for computing the diagonal elements. Our probabilistic bounds for normwise absolute and relative errors apply to Monte Carlo estimators based on random Rademacher, sparse Rademacher, and normalized and unnormalized Gaussian vectors and to vectors with bounded fourth moments. The novel use of matrix concentration inequalities in our proofs represents a systematic model for future analyses. Our bounds mostly do not depend explicitly on the matrix dimension, target different error measures than existing work, and imply that the accuracy of the estimators increases with the diagonal dominance of the matrix. Applications to derivative-based global sensitivity metrics and node centrality measures in network science corroborate this, as do numerical experiments on synthetic test matrices. We recommend against the use in practice of sparse Rademacher vectors, which are the basis for many randomized sketching and sampling algorithms, because they tend to deliver barely a digit of accuracy even under large sampling amounts.}, number={1}, journal={SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS}, author={Hallman, Eric and Ipsen, Ilse C. F. and Saibaba, Arvind K.}, year={2023}, pages={240–269} }
 @article{al daas_ballard_cazeaux_hallman_miedlar_pasha_reid_saibaba_2023, title={RANDOMIZED ALGORITHMS FOR ROUNDING IN THE TENSOR-TRAIN FORMAT}, volume={45}, ISSN={["1095-7197"]}, DOI={10.1137/21M1451191}, abstractNote={The Tensor-Train (TT) format is a highly compact low-rank representation for high-dimensional tensors. TT is particularly useful when representing approximations to the solutions of certain types of parametrized partial differential equations. For many of these problems, computing the solution explicitly would require an infeasible amount of memory and computational time. While the TT format makes these problems tractable, iterative techniques for solving the PDEs must be adapted to perform arithmetic while maintaining the implicit structure. The fundamental operation used to maintain feasible memory and computational time is called rounding, which truncates the internal ranks of a tensor already in TT format. We propose several randomized algorithms for this task that are generalizations of randomized low-rank matrix approximation algorithms and provide significant reduction in computation compared to deterministic TT-rounding algorithms. Randomization is particularly effective in the case of rounding a sum of TT-tensors (where we observe 20x speedup), which is the bottleneck computation in the adaptation of GMRES to vectors in TT format. We present the randomized algorithms and compare their empirical accuracy and computational time with deterministic alternatives.}, number={1}, journal={SIAM JOURNAL ON SCIENTIFIC COMPUTING}, author={Al Daas, Hussam and Ballard, Grey and Cazeaux, Paul and Hallman, Eric and Miedlar, Agnieszka and Pasha, Mirjeta and Reid, Tim W. and Saibaba, Arvind K.}, year={2023}, pages={A74–A95} }
 @article{antil_saibaba_2023, title={Randomized reduced basis methods for parameterized fractional elliptic PDEs}, volume={227}, ISSN={["1872-6925"]}, DOI={10.1016/j.finel.2023.104046}, abstractNote={This paper is interested in developing reduced order models (ROMs) for repeated simulation of fractional elliptic partial differential equations (PDEs) for multiple values of the parameters (e.g., diffusion coefficients or fractional exponent) governing these models. These problems arise in many applications including simulating Gaussian processes, geophysical electromagnetics. The approach uses the Kato integral formula to express the solution as an integral involving the solution of a parameterized elliptic PDE, which is discretized using finite elements in space and sinc quadrature for the fractional part. The offline stage of the ROM is accelerated using a solver for shifted linear systems, MPGMRES-Sh, and using a randomized approach for compressing the snapshot matrix. Our approach is both computational and memory efficient. Numerical experiments on a range of model problems, including an application to Gaussian processes, show the benefits of our approach.}, journal={FINITE ELEMENTS IN ANALYSIS AND DESIGN}, author={Antil, Harbir and Saibaba, Arvind K.}, year={2023}, month={Dec} }
 @article{farazmand_saibaba_2023, title={Tensor-based flow reconstruction from optimally located sensor measurements}, volume={962}, ISSN={["1469-7645"]}, DOI={10.1017/jfm.2023.269}, abstractNote={Reconstructing high-resolution flow fields from sparse measurements is a major challenge in fluid dynamics. Existing methods often vectorize the flow by stacking different spatial directions on top of each other, hence confounding the information encoded in different dimensions. Here, we introduce a tensor-based sensor placement and flow reconstruction method which retains and exploits the inherent multidimensionality of the flow. We derive estimates for the flow reconstruction error, storage requirements and computational cost of our method. We show, with examples, that our tensor-based method is significantly more accurate than similar vectorized methods. Furthermore, the variance of the error is smaller when using our tensor-based method. While the computational cost of our method is comparable to similar vectorized methods, it reduces the storage cost by several orders of magnitude. The reduced storage cost becomes even more pronounced as the dimension of the flow increases. We demonstrate the efficacy of our method on three examples: a chaotic Kolmogorov flow, in situ and satellite measurements of the global sea surface temperature and three-dimensional unsteady simulated flow around a marine research vessel.}, journal={JOURNAL OF FLUID MECHANICS}, author={Farazmand, Mohammad and Saibaba, Arvind K.}, year={2023}, month={May} }
 @article{cho_chung_miller_saibaba_2022, title={Computationally efficient methods for large-scale atmospheric inverse modeling}, volume={15}, ISSN={["1991-9603"]}, DOI={10.5194/gmd-15-5547-2022}, abstractNote={Abstract. Atmospheric inverse modeling describes the process of estimating greenhouse gas fluxes or air pollution emissions at the Earth's surface using observations of these gases collected in the atmosphere. The launch of new satellites, the expansion of surface observation networks, and a desire for more detailed maps of surface fluxes have yielded numerous computational and statistical challenges for standard inverse modeling frameworks that were often originally designed with much smaller data sets in mind. In this article, we discuss computationally efficient methods for large-scale atmospheric inverse modeling and focus on addressing some of the main computational and practical challenges. We develop generalized hybrid projection methods, which are iterative methods for solving large-scale inverse problems, and specifically we focus on the case of estimating surface fluxes. These algorithms confer several advantages. They are efficient, in part because they converge quickly, they exploit efficient matrix–vector multiplications, and they do not require inversion of any matrices. These methods are also robust because they can accurately reconstruct surface fluxes, they are automatic since regularization or covariance matrix parameters and stopping criteria can be determined as part of the iterative algorithm, and they are flexible because they can be paired with many different types of atmospheric models. We demonstrate the benefits of generalized hybrid methods with a case study from NASA's Orbiting Carbon Observatory 2 (OCO-2) satellite. We then address the more challenging problem of solving the inverse model when the mean of the surface fluxes is not known a priori; we do so by reformulating the problem, thereby extending the applicability of hybrid projection methods to include hierarchical priors. We further show that by exploiting mathematical relations provided by the generalized hybrid method, we can efficiently calculate an approximate posterior variance, thereby providing uncertainty information.
                    }, number={14}, journal={GEOSCIENTIFIC MODEL DEVELOPMENT}, author={Cho, Taewon and Chung, Julianne and Miller, Scot M. and Saibaba, Arvind K.}, year={2022}, month={Jul}, pages={5547–5565} }
 @article{saibaba_minster_kilmer_2022, title={Efficient randomized tensor-based algorithms for function approximation and low-rank kernel interactions}, volume={48}, ISSN={["1572-9044"]}, DOI={10.1007/s10444-022-09979-7}, abstractNote={In this paper, we introduce a method for multivariate function approximation using function evaluations, Chebyshev polynomials, and tensor-based compression techniques via the Tucker format. We develop novel randomized techniques to accomplish the tensor compression, provide a detailed analysis of the computational costs, provide insight into the error of the resulting approximations, and discuss the benefits of the proposed approaches. We also apply the tensor-based function approximation to develop low-rank matrix approximations to kernel matrices that describe pairwise interactions between two sets of points; the resulting low-rank approximations are efficient to compute and store (the complexity is linear in the number of points). We present an adaptive version of the function and kernel approximation that determines an approximation that satisfies a user-specified relative error over a set of random points. We extend our approach to the case where the kernel requires repeated evaluations for many values of (hyper)parameters that govern the kernel. We give detailed numerical experiments on example problems involving multivariate function approximation, low-rank matrix approximations of kernel matrices involving well-separated clusters of sources and target points, and a global low-rank approximation of kernel matrices with an application to Gaussian processes. We observe speedups up to 18X over standard matrix-based approaches.}, number={5}, journal={ADVANCES IN COMPUTATIONAL MATHEMATICS}, author={Saibaba, Arvind K. and Minster, Rachel and Kilmer, Misha E.}, year={2022}, month={Oct} }
 @article{majumder_guan_reich_saibaba_2022, title={Kryging: geostatistical analysis of large-scale datasets using Krylov subspace methods}, volume={32}, ISSN={["1573-1375"]}, DOI={10.1007/s11222-022-10104-3}, abstractNote={Analyzing massive spatial datasets using a Gaussian process model poses computational challenges. This is a problem prevailing heavily in applications such as environmental modeling, ecology, forestry and environmental health. We present a novel approximate inference methodology that uses profile likelihood and Krylov subspace methods to estimate the spatial covariance parameters and makes spatial predictions with uncertainty quantification for point-referenced spatial data. “Kryging” combines Kriging and Krylov subspace methods and applies for both observations on regular grid and irregularly spaced observations, and for any Gaussian process with a stationary isotropic (and certain geometrically anisotropic) covariance function, including the popular Matérn covariance family. We make use of the block Toeplitz structure with Toeplitz blocks of the covariance matrix and use fast Fourier transform methods to bypass the computational and memory bottlenecks of approximating log-determinant and matrix-vector products. We perform extensive simulation studies to show the effectiveness of our model by varying sample sizes, spatial parameter values and sampling designs. A real data application is also performed on a dataset consisting of land surface temperature readings taken by the MODIS satellite. Compared to existing methods, the proposed method performs satisfactorily with much less computation time and better scalability.}, number={5}, journal={STATISTICS AND COMPUTING}, author={Majumder, Suman and Guan, Yawen and Reich, Brian J. and Saibaba, Arvind K.}, year={2022}, month={Oct} }
 @article{dudley_saibaba_alexanderian_2022, title={MONTE CARLO ESTIMATORS FOR THE SCHATTEN p-NORM OF SYMMETRIC POSITIVE SEMIDEFINITE MATRICES}, volume={55}, ISSN={["1068-9613"]}, url={http://dx.doi.org/10.1553/etna_vol55s213}, DOI={10.1553/etna_vol55s213}, abstractNote={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.}, journal={ELECTRONIC TRANSACTIONS ON NUMERICAL ANALYSIS}, publisher={Osterreichische Akademie der Wissenschaften, Verlag}, author={DUDLEY, E. T. H. A. N. and SAIBABA, A. R. V. I. N. D. K. and ALEXANDERIAN, A. L. E. N.}, year={2022}, pages={213–241} }
 @article{kilmer_saibaba_2022, title={STRUCTURED MATRIX APPROXIMATIONS VIA TENSOR DECOMPOSITIONS}, volume={43}, ISSN={["1095-7162"]}, DOI={10.1137/21M1418290}, abstractNote={We provide a computational framework for approximating a class of structured matrices; here, the term structure is very general, and may refer to a regular sparsity pattern (e.g., block-banded), or be more highly structured (e.g., symmetric block Toeplitz). The goal is to uncover {\it additional latent structure} that will in turn lead to computationally efficient algorithms when the new structured matrix approximations are employed in the place of the original operator. Our approach has three steps: map the structured matrix to tensors, use tensor compression algorithms, and map the compressed tensors back to obtain two different matrix representations -- sum of Kronecker products and block low-rank format. The use of tensor decompositions enables us to uncover latent structure in the problem and leads to compressed representations of the original matrix that can be used efficiently in applications. The resulting matrix approximations are memory efficient, easy to compute with, and preserve the error that is due to the tensor compression in the Frobenius norm. Our framework is quite general. We illustrate the ability of our method to uncover block-low-rank format on structured matrices from two applications: system identification, space-time covariance matrices. In addition, we demonstrate that our approach can uncover sum of structured Kronecker products structure on several matrices from the SuiteSparse collection. Finally, we show that our framework is broad enough to encompass and improve on other related results from the literature, as we illustrate with the approximation of a three-dimensional blurring operator.}, number={4}, journal={SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS}, author={Kilmer, Misha E. and Saibaba, Arvind K.}, year={2022}, pages={1599–1626} }
 @article{minster_saibaba_kar_chakrabortty_2021, title={EFFICIENT ALGORITHMS FOR EIGENSYSTEM REALIZATION USING RANDOMIZED SVD}, volume={42}, ISSN={["1095-7162"]}, DOI={10.1137/20M1327616}, abstractNote={Eigensystem Realization Algorithm (ERA) is a data-driven approach for subspace system identification and is widely used in many areas of engineering. However, the computational cost of the ERA is dominated by a step that involves the singular value decomposition (SVD) of a large, dense matrix with block Hankel structure. This paper develops computationally efficient algorithms for reducing the computational cost of the SVD step by using randomized subspace iteration and exploiting the block Hankel structure of the matrix. We provide a detailed analysis of the error in the identified system matrices and the computational cost of the proposed algorithms. We demonstrate the accuracy and computational benefits of our algorithms on two test problems: the first involves a partial differential equation that models the cooling of steel rails, and the second is an application from power systems engineering.}, number={2}, journal={SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS}, author={Minster, Rachel and Saibaba, Arvind K. and Kar, Jishnudeep and Chakrabortty, Aranya}, year={2021}, pages={1045–1072} }
 @article{saibaba_hart_bloemen waanders_2021, title={Randomized algorithms for generalized singular value decomposition with application to sensitivity analysis}, volume={28}, ISSN={["1099-1506"]}, DOI={10.1002/nla.2364}, abstractNote={Abstract}, number={4}, journal={NUMERICAL LINEAR ALGEBRA WITH APPLICATIONS}, author={Saibaba, Arvind K. and Hart, Joseph and Bloemen Waanders, Bart}, year={2021}, month={Aug} }
 @article{saibaba_prasad_sturler_miller_kilmer_2021, title={Randomized approaches to accelerate MCMC algorithms for Bayesian inverse problems}, volume={440}, ISSN={["1090-2716"]}, DOI={10.1016/j.jcp.2021.110391}, abstractNote={Markov chain Monte Carlo (MCMC) approaches are traditionally used for uncertainty quantification in inverse problems where the physics of the underlying sensor modality is described by a partial differential equation (PDE). However, the use of MCMC algorithms is prohibitively expensive in applications where each log-likelihood evaluation may require hundreds to thousands of PDE solves corresponding to multiple sensors; i.e., spatially distributed sources and receivers perhaps operating at different frequencies or wavelengths depending on the precise application. We show how to mitigate the computational cost of each log-likelihood evaluation by using several randomized techniques and embed these randomized approximations within MCMC algorithms. The resulting MCMC algorithms are computationally efficient methods for quantifying the uncertainty associated with the reconstructed parameters. We demonstrate the accuracy and computational benefits of our proposed algorithms on a model application from diffuse optical tomography where we invert for the spatial distribution of optical absorption.}, journal={JOURNAL OF COMPUTATIONAL PHYSICS}, author={Saibaba, Arvind K. and Prasad, Pranjal and Sturler, Eric and Miller, Eric and Kilmer, Misha E.}, year={2021}, month={Sep} }
 @article{saibaba_chung_petroske_2020, title={Efficient Krylov subspace methods for uncertainty quantification in large Bayesian linear inverse problems}, volume={27}, ISSN={["1099-1506"]}, DOI={10.1002/nla.2325}, abstractNote={Summary}, number={5}, journal={NUMERICAL LINEAR ALGEBRA WITH APPLICATIONS}, author={Saibaba, Arvind K. and Chung, Julianne and Petroske, Katrina}, year={2020}, month={Oct} }
 @article{miller_saibaba_trudeau_mountain_andrews_2020, title={Geostatistical inverse modeling with very large datasets: an example from the Orbiting Carbon Observatory 2 (OCO-2) satellite}, volume={13}, ISSN={["1991-9603"]}, DOI={10.5194/gmd-13-1771-2020}, abstractNote={Abstract. Geostatistical inverse modeling (GIM) has become a common approach to estimating greenhouse gas fluxes at the Earth's surface using atmospheric observations. GIMs are unique relative to other commonly used approaches because they do not require a single emissions inventory or a bottom–up model to serve as an initial guess of the fluxes. Instead, a modeler can incorporate a wide range of environmental, economic, and/or land use data to estimate the fluxes. Traditionally, GIMs have been paired with in situ observations that number in the thousands or tens of thousands. However, the number of available atmospheric greenhouse gas observations has been increasing enormously as the number of satellites, airborne measurement campaigns, and in situ monitoring stations continues to increase. This era of prolific greenhouse gas observations presents computational and statistical challenges for inverse modeling frameworks that have traditionally been paired with a limited number of in situ monitoring sites. In this article, we discuss the challenges of estimating greenhouse gas fluxes using large atmospheric datasets with a particular focus on GIMs. We subsequently discuss several strategies for estimating the fluxes and quantifying uncertainties, strategies that are adapted from hydrology, applied math, or other academic fields and are compatible with a wide variety of atmospheric models. We further evaluate the accuracy and computational burden of each strategy using a synthetic CO2 case study based upon NASA's Orbiting Carbon Observatory 2 (OCO-2) satellite. Specifically, we simultaneously estimate a full year of 3-hourly CO2 fluxes across North America in one case study – a total of 9.4×106 unknown fluxes using 9.9×104 synthetic observations. The strategies discussed here provide accurate estimates of CO2 fluxes that are comparable to fluxes calculated directly or analytically. We are also able to approximate posterior uncertainties in the fluxes, but these approximations are, typically, an over- or underestimate depending upon the strategy employed and the degree of approximation required to make the calculations manageable.}, number={3}, journal={GEOSCIENTIFIC MODEL DEVELOPMENT}, author={Miller, Scot M. and Saibaba, Arvind K. and Trudeau, Michael E. and Mountain, Marikate E. and Andrews, Arlyn E.}, year={2020}, month={Apr}, pages={1771–1785} }
 @article{herman_alexanderian_saibaba_2020, title={RANDOMIZATION AND REWEIGHTED l(1)-MINIMIZATION FOR A-OPTIMAL DESIGN OF LINEAR INVERSE PROBLEMS}, volume={42}, ISSN={["1095-7197"]}, DOI={10.1137/19M1267362}, abstractNote={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 ...}, number={3}, journal={SIAM JOURNAL ON SCIENTIFIC COMPUTING}, author={Herman, Elizabeth and Alexanderian, Alen and Saibaba, Arvind K.}, year={2020}, pages={A1714–A1740} }
 @article{minster_saibaba_kilmer_2020, title={Randomized Algorithms for Low-Rank Tensor Decompositions in the Tucker Format}, volume={2}, ISSN={["2577-0187"]}, url={http://dx.doi.org/10.1137/19m1261043}, DOI={10.1137/19m1261043}, abstractNote={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.}, number={1}, journal={SIAM JOURNAL ON MATHEMATICS OF DATA SCIENCE}, author={Minster, Rachel and Saibaba, Arvind K. and Kilmer, Misha E.}, year={2020}, pages={189–215} }
 @article{randomized discrete empirical interpolation method for nonlinear model reduction_2020, url={http://dx.doi.org/10.1137/19m1243270}, DOI={10.1137/19m1243270}, abstractNote={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.}, journal={SIAM Journal on Scientific Computing}, year={2020}, month={Jan} }
 @article{saibaba_bardsley_brown_alexanderian_2019, title={Efficient Marginalization-Based MCMC Methods for Hierarchical Bayesian Inverse Problems}, volume={7}, ISSN={["2166-2525"]}, DOI={10.1137/18M1220625}, abstractNote={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.}, number={3}, journal={SIAM-ASA JOURNAL ON UNCERTAINTY QUANTIFICATION}, author={Saibaba, Arvind K. and Bardsley, Johnathan and Brown, D. Andrew and Alexanderian, Alen}, year={2019}, pages={1105–1131} }
 @article{chi_hu_saibaba_rao_2019, title={Going Off the Grid: Iterative Model Selection for Biclustered Matrix Completion}, volume={28}, ISSN={["1537-2715"]}, url={https://doi.org/10.1080/10618600.2018.1482763}, DOI={10.1080/10618600.2018.1482763}, abstractNote={ABSTRACT 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 to select an appropriate amount of row and column smoothing, namely, to 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, combined with a stochastic Quasi-Newton approach. Supplementary material for this article is available online.}, number={1}, journal={JOURNAL OF COMPUTATIONAL AND GRAPHICAL STATISTICS}, publisher={Informa UK Limited}, author={Chi, Eric C. and Hu, Liuyi and Saibaba, Arvind K. and Rao, Arvind U. K.}, year={2019}, month={Jan}, pages={36–47} }
 @article{saibaba_2019, title={RANDOMIZED SUBSPACE ITERATION: ANALYSIS OF CANONICAL ANGLES AND UNITARILY INVARIANT NORMS}, volume={40}, ISSN={["1095-7162"]}, DOI={10.1137/18M1179432}, abstractNote={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.}, number={1}, journal={SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS}, author={Saibaba, Arvind K.}, year={2019}, pages={23–48} }
 @article{zhang_saibaba_kilmer_aeron_2018, title={A randomized tensor singular value decomposition based on the t-product}, volume={25}, ISSN={["1099-1506"]}, DOI={10.1002/nla.2179}, abstractNote={Summary}, number={5}, journal={NUMERICAL LINEAR ALGEBRA WITH APPLICATIONS}, author={Zhang, Jiani and Saibaba, Arvind K. and Kilmer, Misha E. and Aeron, Shuchin}, year={2018}, month={Oct} }
 @article{alexanderian_saibaba_2018, title={EFFICIENT D-OPTIMAL DESIGN OF EXPERIMENTS FOR INFINITE-DIMENSIONAL BAYESIAN LINEAR INVERSE PROBLEMS}, volume={40}, ISSN={["1095-7197"]}, DOI={10.1137/17M115712X}, abstractNote={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.}, number={5}, journal={SIAM JOURNAL ON SCIENTIFIC COMPUTING}, author={Alexanderian, Alen and Saibaba, Arvind K.}, year={2018}, pages={A2956–A2985} }
 @article{chung_saibaba_brown_westman_2018, title={Efficient generalized Golub-Kahan based methods for dynamic inverse problems}, volume={34}, url={http://www.scopus.com/inward/record.url?eid=2-s2.0-85041011795&partnerID=MN8TOARS}, DOI={10.1088/1361-6420/aaa0e1}, abstractNote={We consider efficient methods for computing solutions to and estimating uncertainties in dynamic inverse problems, where the parameters of interest may change during the measurement procedure. Compared to static inverse problems, incorporating prior information in both space and time in a Bayesian framework can become computationally intensive, in part, due to the large number of unknown parameters. In these problems, explicit computation of the square root and/or inverse of the prior covariance matrix is not possible, so we consider efficient, iterative, matrix-free methods based on the generalized Golub–Kahan bidiagonalization that allow automatic regularization parameter and variance estimation. We demonstrate that these methods for dynamic inversion can be more flexible than standard methods and develop efficient implementations that can exploit structure in the prior, as well as possible structure in the forward model. Numerical examples from photoacoustic tomography, space-time deblurring, and passive seismic tomography demonstrate the range of applicability and effectiveness of the described approaches. Specifically, in passive seismic tomography, we demonstrate our approach on both synthetic and real data. To demonstrate the scalability of our algorithm, we solve a dynamic inverse problem with approximately 43000 measurements and 7.8 million unknowns in under 40 s on a standard desktop.}, number={2}, journal={Inverse Problems}, author={Chung, J. and Saibaba, Arvind and Brown, M. and Westman, E.}, year={2018} }
 @article{attia_alexanderian_saibaba_2018, title={Goal-oriented optimal design of experiments for large-scale Bayesian linear inverse problems}, volume={34}, ISSN={["1361-6420"]}, DOI={10.1088/1361-6420/aad210}, abstractNote={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 -norm penalty approach, and propose a practical strategy for specifying the associated penalty parameter.}, number={9}, journal={INVERSE PROBLEMS}, author={Attia, Ahmed and Alexanderian, Alen and Saibaba, Arvind K.}, year={2018}, month={Sep} }
 @article{brown_saibaba_vallelian_2018, title={Low-Rank Independence Samplers in Hierarchical Bayesian Inverse Problems}, volume={6}, ISSN={["2166-2525"]}, DOI={10.1137/17M1137218}, abstractNote={In Bayesian inverse problems, the posterior distribution is used to quantify uncertainty about the reconstructed solution. In fully Bayesian approaches in which prior parameters are assigned hyperpriors, 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.}, number={3}, journal={SIAM-ASA JOURNAL ON UNCERTAINTY QUANTIFICATION}, author={Brown, D. Andrew and Saibaba, Arvind and Vallelian, Sarah}, year={2018}, pages={1076–1100} }
 @article{drmac_saibaba_2018, title={THE DISCRETE EMPIRICAL INTERPOLATION METHOD: CANONICAL STRUCTURE AND FORMULATION IN WEIGHTED INNER PRODUCT SPACES}, volume={39}, ISSN={["1095-7162"]}, DOI={10.1137/17M1129635}, abstractNote={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.}, number={3}, journal={SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS}, author={Drmac, Zlatko and Saibaba, Arvind Krishna}, year={2018}, pages={1152–1180} }
 @article{chung_saibaba_2017, title={GENERALIZED HYBRID ITERATIVE METHODS FOR LARGE-SCALE BAYESIAN INVERSE PROBLEMS}, volume={39}, ISSN={["1095-7197"]}, DOI={10.1137/16m1081968}, abstractNote={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.}, number={5}, journal={SIAM JOURNAL ON SCIENTIFIC COMPUTING}, publisher={Society for Industrial & Applied Mathematics (SIAM)}, author={Chung, Julianne and Saibaba, Arvind K.}, year={2017}, pages={S24–S46} }
 @article{bakhos_kitanidis_ladenheim_saibaba_szyld_2017, title={MULTIPRECONDITIONED GMRES FOR SHIFTED SYSTEMS}, volume={39}, ISSN={["1095-7197"]}, DOI={10.1137/16m1068694}, abstractNote={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.}, number={5}, journal={SIAM JOURNAL ON SCIENTIFIC COMPUTING}, publisher={Society for Industrial & Applied Mathematics (SIAM)}, author={Bakhos, Tania and Kitanidis, Peter K. and Ladenheim, Scott and Saibaba, Arvind K. and Szyld, Daniel B.}, year={2017}, pages={S222–S247} }
 @article{saibaba_alexanderian_ipsen_2017, title={Randomized matrix-free trace and log-determinant estimators}, volume={137}, ISSN={["0945-3245"]}, url={http://www.scopus.com/inward/record.url?eid=2-s2.0-85017115710&partnerID=MN8TOARS}, DOI={10.1007/s00211-017-0880-z}, abstractNote={We present randomized algorithms for estimating the trace and determinant 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.}, number={2}, journal={NUMERISCHE MATHEMATIK}, author={Saibaba, Arvind K. and Alexanderian, Alen and Ipsen, Ilse C. F.}, year={2017}, month={Oct}, pages={353–395} }
 @article{saibaba_2016, title={HOID: HIGHER ORDER INTERPOLATORY DECOMPOSITION FOR TENSORS BASED ON TUCKER REPRESENTATION}, volume={37}, ISSN={["1095-7162"]}, url={http://www.scopus.com/inward/record.url?eid=2-s2.0-84990853624&partnerID=MN8TOARS}, DOI={10.1137/15m1048628}, abstractNote={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.}, number={3}, journal={SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS}, author={Saibaba, Arvind K.}, year={2016}, pages={1223–1249} }
 @inproceedings{saibaba_krishnamurthy_anderson_kainerstorfer_sassaroli_miller_fantini_kilmer_2015, title={3D parameter reconstruction in hyperspectral diffuse optical tomography}, volume={9319}, url={http://www.scopus.com/inward/record.url?eid=2-s2.0-84928719400&partnerID=MN8TOARS}, DOI={10.1117/12.2079775}, abstractNote={The imaging of shape perturbation and chromophore concentration using Diffuse Optical Tomography (DOT) data can be mathematically described as an ill-posed and non-linear inverse problem. The reconstruction algorithm for hyperspectral data using a linearized Born model is prohibitively expensive, both in terms of computation and memory. We model the shape of the perturbation using parametric level-set approach (PaLS). We discuss novel computational strategies for reducing the computational cost based on a Krylov subspace approach for parameteric linear systems and a compression strategy for the parameter-to-observation map. We will demonstrate the validity of our approach by comparison with experiments.}, booktitle={Progress in Biomedical Optics and Imaging - Proceedings of SPIE}, author={Saibaba, A.K. and Krishnamurthy, N. and Anderson, P.G. and Kainerstorfer, J.M. and Sassaroli, A. and Miller, E.L. and Fantini, S. and Kilmer, M.E.}, year={2015} }
 @article{bakhos_saibaba_kitanidis_2015, title={A fast algorithm for parabolic PDE-based inverse problems based on Laplace transforms and flexible Krylov solvers}, volume={299}, url={http://www.scopus.com/inward/record.url?eid=2-s2.0-84938709974&partnerID=MN8TOARS}, DOI={10.1016/j.jcp.2015.07.007}, abstractNote={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.}, journal={Journal of Computational Physics}, author={Bakhos, T. and Saibaba, A.K. and Kitanidis, P.K.}, year={2015}, pages={940–954} }
 @article{saibaba_miller_kitanidis_2015, title={Fast Kalman filter using hierarchical matrices and a low-rank perturbative approach}, volume={31}, url={http://www.scopus.com/inward/record.url?eid=2-s2.0-84920486586&partnerID=MN8TOARS}, DOI={10.1088/0266-5611/31/1/015009}, abstractNote={We develop a fast algorithm for a 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 preceding terms and the computation of eigenmodes of the generalized eigenvalue problem. The resulting algorithm for the Kalman filter with a random walk forecast model scales as  ( N ) ?> in memory and  ( 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 CO2 in the subsurface using travel-time tomography.}, number={1}, journal={Inverse Problems}, author={Saibaba, A.K. and Miller, E.L. and Kitanidis, P.K.}, year={2015} }
 @article{saibaba_kilmer_miller_fantini_2015, title={Fast algorithms for hyperspectral diffuse optical tomography}, volume={37}, url={http://www.scopus.com/inward/record.url?eid=2-s2.0-84945917007&partnerID=MN8TOARS}, DOI={10.1137/140990073}, abstractNote={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.}, number={5}, journal={SIAM Journal on Scientific Computing}, author={Saibaba, A.K. and Kilmer, M. and Miller, E.L. and Fantini, S.}, year={2015}, pages={B712–B743} }
 @article{saibaba_kitanidis_2015, title={Fast computation of uncertainty quantification measures in the geostatistical approach to solve inverse problems}, volume={82}, url={http://www.scopus.com/inward/record.url?eid=2-s2.0-84929572878&partnerID=MN8TOARS}, DOI={10.1016/j.advwatres.2015.04.012}, abstractNote={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.}, journal={Advances in Water Resources}, author={Saibaba, A.K. and Kitanidis, P.K.}, year={2015}, pages={124–138} }
 @article{saibaba_lee_kitanidis_2015, title={Randomized algorithms for generalized Hermitian eigenvalue problems with application to computing Karhunen-Loève expansion}, volume={23}, ISSN={1070-5325}, url={http://dx.doi.org/10.1002/nla.2026}, DOI={10.1002/nla.2026}, abstractNote={Summary}, number={2}, journal={Numerical Linear Algebra with Applications}, publisher={Wiley}, author={Saibaba, Arvind K. and Lee, Jonghyun and Kitanidis, Peter K.}, year={2015}, month={Nov}, pages={314–339} }
 @inproceedings{saibaba_miller_kitandis_2014, title={A fast Kalman filter for time-lapse electrical resistivity tomography}, url={http://www.scopus.com/inward/record.url?eid=2-s2.0-84911398513&partnerID=MN8TOARS}, DOI={10.1109/IGARSS.2014.6947146}, abstractNote={We present a reduced complexity algorithm for time-lapse Electrical Resistivity Tomography (ERT) based on an extended Kalman filter. The key idea of the fast algorithm is an efficient representation of state covariance matrix at each step as a weighted combination of the system noise covariance matrix and a low-rank perturbation term. We propose an efficient algorithm for updating the weights and the basis of the low-rank perturbation. The overall computational cost at each iteration is O(Nnm) and storage cost O(N), where N is the number of grid points, and nm is the number of measurements. The performance of this algorithm is demonstrated on a challenging application of monitoring the CO2 plume using synthetic ERT data.}, booktitle={International Geoscience and Remote Sensing Symposium (IGARSS)}, author={Saibaba, A.K. and Miller, E.L. and Kitandis, P.K.}, year={2014}, pages={3152–3155} }
 @article{saibaba_bakhos_kitanidis_2013, title={A flexible krylov solver for shifted systemswith application to oscillatory hydraulic tomography}, volume={35}, url={http://www.scopus.com/inward/record.url?eid=2-s2.0-84892607644&partnerID=MN8TOARS}, DOI={10.1137/120902690}, abstractNote={We discuss efficient solutions to systems of shifted linear systems arising in computations for oscillatory hydraulic tomography. 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 [Kitanidis, Water Resources Res., 31 (1995), pp. 2411--2419]. 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 [Gu, Zhou, and Lin, J. Comput. Math., 25 (2007), pp. 522--530]. 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 appl...}, number={6}, journal={SIAM Journal on Scientific Computing}, author={Saibaba, A.K. and Bakhos, T. and Kitanidis, P.K.}, year={2013} }
 @article{ambikasaran_saibaba_darve_kitanidis_2013, title={Fast Algorithms for Bayesian Inversion}, DOI={10.1007/978-1-4614-7434-0_5}, abstractNote={In this article, we review a few fast algorithms for solving large-scale stochastic inverse problems using Bayesian methods. After a brief introduction to the Bayesian stochastic inverse methodology, we review the following computational techniques, to solve large scale problems: the fast Fourier transform, the fast multipole method (classical and a black-box version), and finallym the hierarchical matrix approach. We emphasize that this is mainly a survey paper presenting a few fast algorithms applicable to large-scale Bayesian inversion techniques, applicable to applications arising from geostatistics. The article is presented at a level accessible to graduate students and computational engineers. Hence, we mainly present the algorithmic ideas and theoretical results.}, journal={Computational Challenges in the Geosciences}, publisher={Springer New York}, author={Ambikasaran, Sivaram and Saibaba, Arvind K. and Darve, Eric F. and Kitanidis, Peter K.}, year={2013}, pages={101–142} }
 @article{saibaba_ambikasaran_yue li_kitanidis_darve_2012, title={Application of hierarchical matrices to linear inverse problems in geostatistics | Application des matrices hiérarchiques aux problèmes d'inversion linéaire en géostatistique}, volume={67}, url={http://www.scopus.com/inward/record.url?eid=2-s2.0-84872037702&partnerID=MN8TOARS}, DOI={10.2516/ogst/2012064}, number={5}, journal={Oil and Gas Science and Technology}, author={Saibaba, A.K. and Ambikasaran, S. and Yue Li, J. and Kitanidis, P.K. and Darve, E.F.}, year={2012}, pages={857–875} }
 @article{saibaba_kitanidis_2012, title={Efficient methods for large-scale linear inversion using a geostatistical approach}, volume={48}, url={http://www.scopus.com/inward/record.url?eid=2-s2.0-84861357098&partnerID=MN8TOARS}, DOI={10.1029/2011WR011778}, abstractNote={In geophysical inverse problems, such as estimating the unknown parameter field from noisy observations of dependent quantities, e.g., hydraulic conductivity from head observations, stochastic Bayesian and geostatistical approaches are frequently used. To obtain best estimates and conditional realizations it is required to perform several matrix‐matrix computations involving the covariance matrix of the discretized field of the parameters. In realistic three‐dimensional fields that are finely discretized, these operations as performed in conventional algorithms become extremely expensive and even prohibitive in terms of memory and computational requirements. Using Hierarchical Matrices, we show how to reduce the complexity of forming approximate matrix‐vector products involving the Covariance matrices in log linear complexity for an arbitrary distribution of points and a wide variety of generalized covariance functions. The resulting system of equations is solved iteratively using a matrix‐free Krylov subspace approach. Furthermore, we show how to generate unconditional realizations using an approximation to the square root of the covariance matrix using Chebyshev matrix polynomials and use the above to generate conditional realizations. We demonstrate the efficiency of our method on a few standard test problems, such as interpolation from noisy observations and contaminant source identification.}, number={5}, journal={Water Resources Research}, author={Saibaba, A.K. and Kitanidis, P.K.}, year={2012} }