Search

We compare the properties of the stable rank and intrinsic dimension of real and complex matrices to those of the classical rank. Basic proofs and examples illustrate that the stable rank does not satisfy any of the fundamental rank properties, while the intrinsic dimension satisfies a few. In particular, the stable rank and intrinsic dimension of a submatrix can exceed those of the original matrix; adding a Hermitian positive semi-definite matrix can lower the intrinsic dimension of the sum; and multiplication by a nonsingular matrix can drastically change the stable rank and the intrinsic dimension. We generalize the concept of stable rank to the p-stable in a Schatten p-norm, thereby unifying the concepts of stable rank and intrinsic dimension: The stable rank is the 2-stable rank, while the intrinsic dimension is the 1-stable rank of a Hermitian positive semi-definite matrix. We derive sum and product inequalities for the pth root of the p-stable rank, and show that it is well-conditioned in the norm-wise absolute sense. The conditioning improves if the matrix and the perturbation are Hermitian positive semi-definite.

Computing low-rank approximations of kernel matrices is an important problem with many applications in scientific computing and data science. We propose methods to efficiently approximate and store low-rank approximations to kernel matrices that depend on certain hyperparameters. The main idea behind our method is to use multivariate Chebyshev function approximation along with the tensor train decomposition of the coefficient tensor. The computations are in two stages: an offline stage, which dominates the computational cost and is parameter-independent, and an online stage, which is inexpensive and instantiated for specific hyperparameters. A variation of this method addresses the case that the kernel matrix is symmetric and positive semi-definite. The resulting algorithms have linear complexity in terms of the sizes of the kernel matrices. We investigate the efficiency and accuracy of our method on parametric kernel matrices induced by various kernels, such as the Mat\'ern kernel, through various numerical experiments. Our methods have speedups up to $200\times$ in the online time compared to other methods with similar complexity and comparable accuracy., Comment: 29 pages, 3 tables, 2 figures

This paper tackles optimal sensor placement for Bayesian linear inverse problems, a popular version of the more general Optimal Experiment Design (OED) problem, using the D-optimality criterion. This is done by establishing connections between sensor placement and Column Subset Selection Problem (CSSP), which is a well-studied problem in Numerical Linear Algebra (NLA). In particular, we use the Golub-Klema-Stewart (GKS) approach which involves computing the truncated Singular Value Decomposition (SVD) followed by a pivoted QR factorization on the right singular vectors. The algorithms are further accelerated by using randomization to compute the low-rank approximation as well as for sampling the indices. The resulting algorithms are robust, computationally efficient, amenable to parallelization, require virtually no parameter tuning, and come with strong theoretical guarantees. One of the proposed algorithms is also adjoint-free which is beneficial in situations, where the adjoint is expensive to evaluate or is not available. Additionally, we develop a method for data completion without solving the inverse problem. Numerical experiments on model inverse problems involving the heat equation and seismic tomography in two spatial dimensions demonstrate the performance of our approaches., Comment: 31 pages, 9 figures

Strong Constraint 4D Variational (SC-4DVAR) is a data assimilation method that is widely used in climate and weather applications. SC-4DVAR involves solving a minimization problem to compute the maximum a posteriori estimate, which we tackle using the Gauss-Newton method. The computation of the descent direction is expensive since it involves the solution of a large-scale and potentially ill-conditioned linear system, solved using the preconditioned conjugate gradient (PCG) method. To address this cost, we efficiently construct scalable preconditioners using three different randomization techniques, which all rely on a certain low-rank structure involving the Gauss-Newton Hessian. The proposed techniques come with theoretical (probabilistic) guarantees on the condition number, and at the same time, are amenable to parallelization. We also develop an adaptive approach to estimate the sketch size and to choose between the reuse or recomputation of the preconditioner. We demonstrate the performance and effectiveness of our methodology on two representative model problems -- the Burgers and barotropic vorticity equation -- showing a drastic reduction in both the number of PCG iterations and the number of Gauss-Newton Hessian products (after including the preconditioner construction cost)., Comment: 28 pages

We study Bayesian methods for large-scale linear inverse problems, focusing on the challenging task of hyperparameter estimation. Typical hierarchical Bayesian formulations that follow a Markov Chain Monte Carlo approach are possible for small problems with very few hyperparameters but are not computationally feasible for problems with a very large number of unknown parameters. In this work, we describe an empirical Bayesian (EB) method to estimate hyperparameters that maximize the marginal posterior, i.e., the probability density of the hyperparameters conditioned on the data, and then we use the estimated values to compute the posterior of the inverse parameters. For problems where the computation of the square root and inverse of prior covariance matrices are not feasible, we describe an approach based on the generalized Golub-Kahan bidiagonalization to approximate the marginal posterior and seek hyperparameters that minimize the approximate marginal posterior. Numerical results from seismic and atmospheric tomography demonstrate the accuracy, robustness, and potential benefits of the proposed approach.

This paper expands the analysis of randomized low-rank approximation beyond the Gaussian distribution to four classes of random matrices: (1) independent sub-Gaussian entries, (2) independent sub-Gaussian columns, (3) independent bounded columns, and (4) independent columns with bounded second moment. Using a novel interpretation of the low-rank approximation error involving sample covariance matrices, we provide insight into the requirements of a \textit{good random matrix} for the purpose of randomized low-rank approximation. Although our bounds involve unspecified absolute constants (a consequence of the underlying non-asymptotic theory of random matrices), they allow for qualitative comparisons across distributions. The analysis offers some details on the minimal number of samples (the number of columns $\ell$ of the random matrix $\boldsymbol\Omega$) and the error in the resulting low-rank approximation. We illustrate our analysis in the context of the randomized subspace iteration method as a representative algorithm for low-rank approximation, however, all the results are broadly applicable to other low-rank approximation techniques. We conclude our discussion with numerical examples using both synthetic and real-world test matrices., Comment: 37 pages, appendix is really supplementary materials

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, and geophysical electromagnetics. The approach uses the Kato integral formula to express the solution as an integral involving the solution of a parametrized 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., Comment: 15 pages, 5 figures, 2 tables

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 3D unsteady simulated flow around a marine research vessel., Comment: 9 figures, 17 pages. A new three-dimensional example is added

We develop hybrid projection methods for computing solutions to large-scale inverse problems, where the solution represents a sum of different stochastic components. Such scenarios arise in many imaging applications (e.g., anomaly detection in atmospheric emissions tomography) where the reconstructed solution can be represented as a combination of two or more components and each component contains different smoothness or stochastic properties. In a deterministic inversion or inverse modeling framework, these assumptions correspond to different regularization terms for each solution in the sum. Although various prior assumptions can be included in our framework, we focus on the scenario where the solution is a sum of a sparse solution and a smooth solution. For computing solution estimates, we develop hybrid projection methods for solution decomposition that are based on a combined flexible and generalized Golub-Kahan processes. This approach integrates techniques from the generalized Golub-Kahan bidiagonalization and the flexible Krylov methods. The benefits of the proposed methods are that the decomposition of the solution can be done iteratively, and the regularization terms and regularization parameters are adaptively chosen at each iteration. Numerical results from photoacoustic tomography and atmospheric inverse modeling demonstrate the potential for these methods to be used for anomaly detection.

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.

This paper tackles efficient methods for Bayesian inverse problems with priors based on Whittle--Mat\'ern Gaussian random fields. The Whittle--Mat\'ern 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\`eve 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., Comment: 21 pages

There are several challenges associated with inverse problems in which we seek to reconstruct a piecewise constant field, and which we model using multiple level sets. Adopting a Bayesian viewpoint, we impose prior distributions on both the level set functions that determine the piecewise constant regions as well as the parameters that determine their magnitudes. We develop a Gauss-Newton approach with a backtracking line search to efficiently compute the maximum a priori (MAP) estimate as a solution to the inverse problem. We use the Gauss-Newton Laplace approximation to construct a Gaussian approximation of the posterior distribution and use preconditioned Krylov subspace methods to sample from the resulting approximation. To visualize the uncertainty associated with the parameter reconstructions we compute the approximate posterior variance using a matrix-free Monte Carlo diagonal estimator, which we develop in this paper. We will demonstrate the benefits of our approach and solvers on synthetic test problems (photoacoustic and hydraulic tomography, respectively a linear and nonlinear inverse problem) as well as an application to X-ray imaging with real data.

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.

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 have 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., Comment: 27 pages, 7 figures, 4 tables

We consider efficient methods for computing solutions to dynamic inverse problems, where both the quantities of interest and the forward operator (measurement process) may change at different time instances but we want to solve for all the images simultaneously. We are interested in large-scale ill-posed problems that are made more challenging by their dynamic nature and, possibly, by the limited amount of available data per measurement step. To remedy these difficulties, we apply 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 to be solved with a generalized Krylov subspace (GKS) method; the regularization parameter can be defined 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., Comment: 30 pages, 10 figure, 3 tables

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., Comment: 27 pages, 6 figures, 2 tables

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

Analyzing massive spatial datasets using Gaussian process model poses computational challenges. This is a problem prevailing heavily in applications such as environmental modeling, ecology, forestry and environmental heath. 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. The proposed method, Kryging, applies for both observations on regular grid and irregularly-spaced observations, and for any Gaussian process with a stationary covariance function, including the popular $\Matern$ covariance family. We make use of the block Toeplitz structure with Toeplitz blocks of the covariance matrix and use fast Fourier transform methods to alleviate the computational and memory bottlenecks. 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.

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

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.

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., Comment: 24 pages, 10 figures, 1 table

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., Comment: 27 Pages; Accepted for publication in SIAM Journal on Scientific Computing

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

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

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., Comment: 27 pages, 8 figures, 5 tables

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

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., Comment: 25 pages, 4 figures. Accepted, SIAM Journal on Matrix Analysis and Applications

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

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

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. In this work, we develop 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 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, 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 $43,000$ measurements and $7.8$ million unknowns in under $40$ seconds on a standard desktop., Comment: 27 pages, 12 figures

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

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., Comment: 42 pages, 7 figures. Supplementary material (https://github.com/echi/IMS/blob/master/BMC_Supplement_JCGS.pdf) and codes available (https://github.com/echi/IMS)

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.

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.

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., Comment: 38 pages. Accepted for publication in Numerische Mathematik

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., Comment: 24 pages, 4 figures, 1 table

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

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.

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.

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., Comment: published in Inverse Problems, 2015 31 015009

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., Comment: Published in Advances in Water Resources

We describe randomized algorithms for computing the dominant eigenmodes of the Generalized Hermitian Eigenvalue Problem (GHEP) $Ax=\lambda 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\`{e}ve expansion, which is a computationally intensive GHEP problem with rapidly decaying eigenvalues., Comment: Second round review at Numerical Linear Algebra with applications

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.