Adaptivity in both space and time has become the norm for solving problems modeled by partial differential equations. The size of the discretized problem makes uniformly refined grids computationally prohibitive. Adaptive refinement of meshes and time steps allows to capture the phenomena of interest while keeping the cost of a simulation tractable on the current hardware. Many fields in science and engineering require the solution of inverse problems where parameters for a given model are estimated based on available measurement information. In contrast to forward (regular) simulations, inverse problems have not extensively benefited from the adaptive solver technology. Previous research in inverse problems has focused mainly on the continuous approach to calculate sensitivities, and has typically employed fixed time and space meshes in the solution process. Inverse problem solvers that make exclusive use of uniform or static meshes avoid complications such as the differentiation of mesh motion equations, or inconsistencies in the sensitivity equations between subdomains with different refinement levels. However, this comes at the cost of low computational efficiency. More efficient computations are possible through judicious use of adaptive mesh refinement, adaptive time steps, and the discrete adjoint method. This dissertation develops a complete framework for fully discrete adjoint sensitivity analysis and inverse problem solutions, in the context of time dependent, adaptive mesh, and adaptive step models. The discrete framework addresses all the necessary ingredients of a state–of–the–art adaptive inverse solution algorithm: adaptive mesh and time step refinement, solution grid transfer operators, a priori and a posteriori error analysis and estimation, and discrete adjoints for sensitivity analysis of flux–limited numerical algorithms. / Ph. D.

Numerous scientific applications have seen the rise of massive inverse problems, where there are too much data to implement an all-at-once strategy to compute a solution. Additionally, tools for regularizing ill-posed inverse problems are infeasible when the problem is too large. This thesis focuses on the development of row-action methods, which can be used to iteratively solve inverse problems when it is not possible to access the entire data-set or forward model simultaneously. We investigate these techniques for linear inverse problems and for separable, nonlinear inverse problems where the objective function is nonlinear in one set of parameters and linear in another set of parameters. For the linear problem, we perform a convergence analysis of these methods, which shows favorable asymptotic and initial convergence properties, as well as a trade-off between convergence rate and precision of iterates that is based on the step-size. These row-action methods can be interpreted as stochastic Newton and stochastic quasi-Newton approaches on a reformulation of the least squares problem, and they can be analyzed as limited memory variants of the recursive least squares algorithm. For ill-posed problems, we introduce sampled regularization parameter selection techniques, which include sampled variants of the discrepancy principle, the unbiased predictive risk estimator, and the generalized cross-validation. We demonstrate the effectiveness of these methods using examples from super-resolution imaging, tomography reconstruction, and image classification. / Doctor of Philosophy / Numerous scientific problems have seen the rise of massive data sets. An example of this is super-resolution, where many low-resolution images are used to construct a high-resolution image, or 3-D medical imaging where a 3-D image of an object of interest with hundreds of millions voxels is reconstructed from x-rays moving through that object. This work focuses on row-action methods that numerically solve these problems by repeatedly using smaller samples of the data to avoid the computational burden of using the entire data set at once. When data sets contain measurement errors, this can cause the solution to get contaminated with noise. While there are methods to handle this issue, when the data set becomes massive, these methods are no longer feasible. This dissertation develops techniques to avoid getting the solution contaminated with noise, even when the data set is immense. The methods developed in this work are applied to numerous scientific applications including super-resolution imaging, tomography, and image classification.

Probing methods for trace estimation of large, sparse matrices has been studied for several decades. In recent years, there has been some work to extend these techniques to instead estimate the diagonal entries of these systems directly. We extend some analysis of trace estimators to their corresponding diagonal estimators, propose a new class of deterministic diagonal estimators which are well-suited to parallel architectures along with heuristic arguments for the design choices in their construction, and conclude with numerical results on diagonal estimation and ordering problems, demonstrating the strengths of our newly-developed methods alongside existing methods. / Master of Science / In the past several decades, as computational resources increase, a recurring problem is that of estimating certain properties very large linear systems (matrices containing real or complex entries). One particularly important quantity is the trace of a matrix, defined as the sum of the entries along its diagonal. In this thesis, we explore a problem that has only recently been studied, in estimating the diagonal entries of a particular matrix explicitly. For these methods to be computationally more efficient than existing methods, and with favorable convergence properties, we require the matrix in question to have a majority of its entries be zero (the matrix is sparse), with the largest-magnitude entries clustered near and on its diagonal, and very large in size. In fact, this thesis focuses on a class of methods called probing methods, which are of particular efficiency when the matrix is not known explicitly, but rather can only be accessed through matrix vector multiplications with arbitrary vectors. Our contribution is new analysis of these diagonal probing methods which extends the heavily-studied trace estimation problem, new applications for which probing methods are a natural choice for diagonal estimation, and a new class of deterministic probing methods which have favorable properties for large parallel computing architectures which are becoming ever-more-necessary as problem sizes continue to increase beyond the scope of single processor architectures.

This research develops a fault diagnosis method for complex systems in the presence of uncertainties and possibility of multiple solutions. Fault diagnosis is a challenging problem because data used in diagnosis contain random errors and often systematic errors as well. Furthermore, fault diagnosis is basically an inverse problem so that it inherits unfavorable characteristics of inverse problems: The existence and uniqueness of an inverse solution are not guaranteed and the solution may be unstable. The weighted least squares method and its variations are traditionally used for solving inverse problems. However, the existing algorithms often fail to identify multiple solutions if they are present. In addition, the existing algorithms are not capable of selecting variables systematically so that they generally use the full model in which may contain unnecessary variables as well as necessary variables. Ignoring this model uncertainty often gives rise to, so called, the smearing effect in solutions, because of which unnecessary variables are overestimated and necessary variables are underestimated. The proposed method solves the inverse problem using Bayesian inference. An engineering system can be parameterized using state variables. The probability of each state variable is inferred from observations made on the system. A bias in an observation is treated as a variable, and the probability of the bias variable is inferred as well. To take the uncertainty of model structure into account, multiple Bayesian models are created with various combinations of the state variables and the bias variables. The results from all models are averaged according to how likely each model is. Gibbs sampling is used for approximating updated probabilities. The method is demonstrated for two applications: the status matching of a turbojet engine and the fault diagnosis of an industrial gas turbine. In the status matching application only physical faults in the components of a turbojet engine are considered whereas in the fault diagnosis application sensor biases are considered as well as physical faults. The proposed method is tested in various faulty conditions using simulated measurements. Results show that the proposed method identifies physical faults and sensor biases simultaneously. It is also demonstrated that multiple solutions can be identified. Overall, there is a clear improvement in ability to identify correct solutions over the full model that contains all state and bias variables.

This thesis is devoted to solving the inverse source problem arising in image reconstruction problems. In general, the solution is non-unique and the problem is severely ill-posed. Therefore, small perturbations, such as the noise in the data, and the modeling error in the forward problem, will cause huge errors in the computations. In practice, the most widely used method to tackle the problem is based on Tikhonov-type regularizations, which minimizes a cost function combining a regularization term and a data fitting term. However, because the two tasks, namely regularization and data fitting, are coupled together in Tikhonov regularization, they are difficult to solve. It happens even if each task can be efficiently solved when they are separate. We propose a method to overcome the major difficulties, namely the non-uniqueness of the solution and noisy data fitting, separately. First we find a particular solution called the orthogonal solution that satisfies the data fitting term. Then we add to it a correction function in the kernel space so that the final solution fulfills the regularization and other physical requirements. The key idea is that the correction function in the kernel has no impact to the data fitting, and the regularization is imposed in a smaller space. Moreover, there is no parameter needed to balance the data fitting and regularization terms. As a case study, we apply the proposed method to Fluorescence Tomography (FT), an emerging imaging technique well known for its ill-posedness and low image resolution in existing reconstruction techniques. We demonstrate by theory and examples that the proposed algorithm can drastically improve the computation speed and the image resolution over existing methods.

Cauchy problems for elliptic partial differential equations arise in many important applications, such as, cardiography, nondestructive testing, heat transfer, sonic boom produced by a maneuvering aerofoil, etc. Elliptic Cauchy problems are typically ill-posed, i.e., there may not be a solution for some Cauchy data, and even if a solution exists uniquely, it may not depend continuously on the Cauchy data. The ill-posedness causes numerical instability and makes the classical numerical methods inappropriate to solve such problems. For Cauchy problems, the research on uniqueness, stability, and efficient numerical methods are of significant interest to mathematicians. The main focus of this thesis is to develop numerical techniques for elliptic Cauchy problems. Elliptic Cauchy problems can be approached as data completion problems, i.e., from over-specified Cauchy data on an accessible part of the boundary, one can try to recover missing data on the inaccessible part of the boundary. Then, the Cauchy problems can be solved by finding a so-lution to a well-posed boundary value problem for which the recovered data constitute a boundary condition on the inaccessible part of the boundary. In this thesis, we use natural linearization approach to transform the linear Cauchy problem into a problem of solving a linear operator equation. We consider this operator in a weaker image space H−1, which differs from the previous works where the image space of the operator is usually considered as L2 . The lower smoothness of the image space will make a problem a bit more ill-posed. But under such settings, we can prove the compactness of the considered operator. At the same time, it allows a relaxation of the assumption concerning noise. The numerical methods that can cope with these ill-posed operator equations are the so called regularization methods. One prominent example of such regularization methods is Tikhonov regularization which is frequently used in practice. Tikhonov regularization can be considered as a least-squares tracking of data with a regularization term. In this thesis we discuss a possibility to improve the reconstruction accuracy of the Tikhonov regularization method by using an iterative modification of Tikhonov regularization. With this iterated Tikhonov regularization the effect of the penalty term fades away as iterations go on. In the application of iterated Tikhonov regularization, we find that for severely ill-posed problems such as elliptic Cauchy problems, discretization has such a powerful influence on the accuracy of the regularized solution that only with some reasonable discretization level, desirable accuracy can be achieved. Thus, regularization by projection method which is commonly known as self-regularization is also considered in this thesis. With this method, the regularization is achieved only by discretization along with an appropriate choice of discretization level. For all regularization methods, the choice of an appropriate regularization parameter is a crucial issue. For this purpose, we propose the balancing principle which is a recently introduced powerful technique for the choice of the regularization parameter. While applying this principle, a balance between the components related to the convergence rate and stability in the accuracy estimates has to be made. The main advantage of the balancing principle is that it can work in an adaptive way to obtain an appropriate value of the regularization parameter, and it does not use any quantitative knowledge of convergence rate or stability. The accuracy provided by this adaptive strategy is worse only by a constant factor than one could achieve in the case of known stability and convergence rates. We apply the balancing principle in both iterated Tikhonov regularization and self-regularization methods to choose the proper regularization parameters. In the thesis, we also investigate numerical techniques based on iterative Tikhonov regular-ization for nonlinear elliptic Cauchy problems. We consider two types of problems. In the first kind, the nonlinear problem can be transformed to a linear problem while in the second kind, linearization of the nonlinear problem is not possible, and for this we propose a special iterative method which differs from methods such as Landweber iteration and Newton-type method which are usually based on the calculation of the Frech´et derivative or adjoint of the equation. Abundant examples are presented in the thesis, which illustrate the performance of the pro-posed regularization methods as well as the balancing principle. At the same time, these examples can be viewed as a support for the theoretical results achieved in this thesis. In the end of this thesis, we describe the sonic boom problem, where we first encountered the ill-posed nonlinear Cauchy problem. This is a very difficult problem and hence we took this problem to provide a motivation for the model problems. These model problems are discussed one by one in the thesis in the increasing order of difficulty, ending with the nonlinear problems in Chapter 5. The main results of the dissertation are communicated in the article [35].

Two inverse problems in elasticity are addressed with motivation from cellular biomechanics. The first application is computation of holding forces on a cell during its manipulation and the second application is estimation of a cell’s interior elastic mapping (i.e., estimation of inhomogeneous distribution of stiffness) using only boundary forces and displacements. It is clear from recent works that mechanical forces can play an important role in developmental biology. In this regard, we have developed a vision-based force-sensing technique to estimate forces that are acting on a cell while it is manipulated. This problem is connected to one inverse problem in elasticity known as Cauchy’s problem in elasticity. Geometric nonlinearity under noisy displacement data is accounted while developing the solution procedures for Cauchy’s problem. We have presented solution procedures to the Cauchy’s problem under noisy displacement data. Geometric nonlinearity is also considered in order to account large deformations that the mechanisms (grippers) undergo during the manipulation. The second inverse problem is connected to elastic mapping of the cell. We note that recent works in biomechanics have shown that the disease state can alter the gross stiffness of a cell. Therefore, the pertinent question that one can ask is that which portion (for example Nucleus, cortex, ER) of the elastic property of the cell is majorly altered by the disease state. Mathematically, this question (estimation of inhomogeneous properties of cell) can be answered by solving an inverse elastic boundary value problem using sets of force-displacements boundary measurements. We address the theoretical question of number of boundary data sets required to solve the inverse boundary value problem.

For many scientific applications, inverse problems have played a key role in solving important problems by enabling researchers to estimate desired parameters of a system from observed measurements. For example, large-scale inverse problems arise in many global problems and medical imaging problems such as greenhouse gas tracking and computational tomography reconstruction. This dissertation describes advancements in computational tools for solving large-scale inverse problems and for uncertainty quantification. Oftentimes, inverse problems are ill-posed and large-scale. Iterative projection methods have dramatically reduced the computational costs of solving large-scale inverse problems, and regularization methods have been critical in obtaining stable estimations by applying prior information of unknowns via Bayesian inference. However, by combining iterative projection methods and variational regularization methods, hybrid projection approaches, in particular generalized hybrid methods, create a powerful framework that can maximize the benefits of each method. In this dissertation, we describe various advancements and extensions of hybrid projection methods that we developed to address three recent open problems. First, we develop hybrid projection methods that incorporate mixed Gaussian priors, where we seek more sophisticated estimations where the unknowns can be treated as random variables from a mixture of distributions. Second, we describe hybrid projection methods for mean estimation in a hierarchical Bayesian approach. By including more than one prior covariance matrix (e.g., mixed Gaussian priors) or estimating unknowns and hyper-parameters simultaneously (e.g., hierarchical Gaussian priors), we show that better estimations can be obtained. Third, we develop computational tools for a respirometry system that incorporate various regularization methods for both linear and nonlinear respirometry inversions. For the nonlinear systems, blind deconvolution methods are developed and prior knowledge of nonlinear parameters are used to reduce the dimension of the nonlinear systems. Simulated and real-data experiments of the respirometry problems are provided. This dissertation provides advanced tools for computational inversion and uncertainty quantification. / Doctor of Philosophy / For many scientific applications, inverse problems have played a key role in solving important problems by enabling researchers to estimate desired parameters of a system from observed measurements. For example, large-scale inverse problems arise in many global problems such as greenhouse gas tracking where the problem of estimating the amount of added or removed greenhouse gas at the atmosphere gets more difficult. The number of observations has been increased with improvements in measurement technologies (e.g., satellite). Therefore, the inverse problems become large-scale and they are computationally hard to solve. Another example of an inverse problem arises in tomography, where the goal is to examine materials deep underground (e.g., to look for gas or oil) or reconstruct an image of the interior of the human body from exterior measurements (e.g., to look for tumors). For tomography applications, there are typically fewer measurements than unknowns, which results in non-unique solutions. In this dissertation, we treat unknowns as random variables with prior probability distributions in order to compensate for a deficiency in measurements. We consider various additional assumptions on the prior distribution and develop efficient and robust numerical methods for solving inverse problems and for performing uncertainty quantification. We apply our developed methods to many numerical applications such as greenhouse gas tracking, seismic tomography, spherical tomography problems, and the estimation of CO2 of living organisms.

Compressed Sensing (CS) is an emerging paradigm in which signals are recovered from undersampled nonadaptive linear measurements taken at a rate proportional to the signal's true information content as opposed to its ambient dimension. The resulting problem consists in finding a sparse solution to an underdetermined system of linear equations. It has now been established, both theoretically and empirically, that certain optimization algorithms are able to solve such problems. Iterative Hard Thresholding (IHT) (Blumensath and Davies, 2007), which is the focus of this thesis, is an established CS recovery algorithm which is known to be effective in practice, both in terms of recovery performance and computational efficiency. However, theoretical analysis of IHT to date suffers from two drawbacks: state-of-the-art worst-case recovery conditions have not yet been quantified in terms of the sparsity/undersampling trade-off, and also there is a need for average-case analysis in order to understand the behaviour of the algorithm in practice. In this thesis, we present a new recovery analysis of IHT, which considers the fixed points of the algorithm. In the context of arbitrary matrices, we derive a condition guaranteeing convergence of IHT to a fixed point, and a condition guaranteeing that all fixed points are 'close' to the underlying signal. If both conditions are satisfied, signal recovery is therefore guaranteed. Next, we analyse these conditions in the case of Gaussian measurement matrices, exploiting the realistic average-case assumption that the underlying signal and measurement matrix are independent. We obtain asymptotic phase transitions in a proportional-dimensional framework, quantifying the sparsity/undersampling trade-off for which recovery is guaranteed. By generalizing the notion of xed points, we extend our analysis to the variable stepsize Normalised IHT (NIHT) (Blumensath and Davies, 2010). For both stepsize schemes, comparison with previous results within this framework shows a substantial quantitative improvement. We also extend our analysis to a related algorithm which exploits the assumption that the underlying signal exhibits tree-structured sparsity in a wavelet basis (Baraniuk et al., 2010). We obtain recovery conditions for Gaussian matrices in a simplified proportional-dimensional asymptotic, deriving bounds on the oversampling rate relative to the sparsity for which recovery is guaranteed. Our results, which are the first in the phase transition framework for tree-based CS, show a further significant improvement over results for the standard sparsity model. We also propose a dynamic programming algorithm which is guaranteed to compute an exact tree projection in low-order polynomial time.

No description available.