Search

The aim of Electrical Impedance Tomography (EIT) is to determine the electrical conductivity distribution inside a domain by applying currents and measuring voltages on its boundary. Mathematically, the EIT reconstruction task can be formulated as a non-linear inverse problem. The Bayesian inverse problems framework has been applied expensively to solutions of the EIT inverse problem, in particular in the cases when the unknown conductivity is believed to be blocky. Recently, the Sparsity Promoting Iterative Alternating Sequential (PS-IAS) algorithm, originally proposed for the solution of linear inverse problems, has been adapted for the non linear case of EIT reconstruction in a computationally efficient manner. Here we introduce a hybrid version of the SP-IAS algorithms for the nonlinear EIT inverse problem, providing a detailed description of the implementation details, with a specific focus on parameters selection. The method is applied to the 2023 Kuopio Tomography Challenge dataset, with a comprehensive report of the running times for the different cases and parameter selections.

The electrical impedance tomography (EIT) problem of estimating the unknown conductivity distribution inside a domain from boundary current or voltage measurements requires the solution of a nonlinear inverse problem. Sparsity promoting hierarchical Bayesian models have been shown to be very effective in the recovery of almost piecewise constant solutions in linear inverse problems. We demonstrate that by exploiting linear algebraic considerations it is possible to organize the calculation for the Bayesian solution of the nonlinear EIT inverse problem via finite element methods with sparsity promoting priors in a computationally efficient manner. The proposed approach uses the Iterative Alternating Sequential (IAS) algorithm for the solution of the linearized problems. Within the IAS algorithm, a substantial reduction in computational complexity is attained by exploiting the low dimensionality of the data space and an adjoint formulation of the Tikhonov regularized solution that constitutes part of the iterative updating scheme. Numerical tests illustrate the computational efficiency of the proposed algorithm. The paper sheds light also on the convexity properties of the objective function of the maximum a posteriori (MAP) estimation problem.

We exploit the similarities between Tikhonov regularization and Bayesian hierarchical models to propose a regularization scheme that acts like a distributed Tikhonov regularization where the amount of regularization varies from component to component. In the standard formulation, Tikhonov regularization compensates for the inherent ill-conditioning of linear inverse problems by augmenting the data fidelity term measuring the mismatch between the data and the model output with a scaled penalty functional. The selection of the scaling is the core problem in Tikhonov regularization. If an estimate of the amount of noise in the data is available, a popular way is to use the Morozov discrepancy principle, stating that the scaling parameter should be chosen so as to guarantee that the norm of the data fitting error is approximately equal to the norm of the noise in the data. A too small value of the regularization parameter would yield a solution that fits to the noise while a too large value would lead to an excessive penalization of the solution. In many applications, it would be preferable to apply distributed regularization, replacing the regularization scalar by a vector valued parameter, allowing different regularization for different components of the unknown, or for groups of them. A distributed Tikhonov-inspired regularization is particularly well suited when the data have significantly different sensitivity to different components, or to promote sparsity of the solution. The numerical scheme that we propose, while exploiting the Bayesian interpretation of the inverse problem and identifying the Tikhonov regularization with the Maximum A Posteriori (MAP) estimation, requires no statistical tools. A combination of numerical linear algebra and optimization tools makes the scheme computationally efficient and suitable for problems where the matrix is not explicitly available.

We consider inverse problems estimating distributed parameters from indirect noisy observations through discretization of continuum models described by partial differential or integral equations. It is well understood that the errors arising from the discretization can be detrimental for ill-posed inverse problems, as discretization error behaves as correlated noise. While this problem can be avoided with a discretization fine enough to suppress the modeling error level below that of the exogenous noise that is addressed, e.g., by regularization, the computational resources needed to deal with the additional degrees of freedom may require high performance computing environment. Following an earlier idea, we advocate the notion that the discretization is one of the unknowns of the inverse problem, and is updated iteratively together with the solution. In this approach, the discretization, defined in terms of an underlying metric, is refined selectively only where the representation power of the current mesh is insufficient. In this paper we allow the metrics and meshes to be anisotropic, and we show that this leads to significant reduction of memory allocation and computing time.

Dictionary learning methods continue to gain popularity for the solution of challenging inverse problems. In the dictionary learning approach, the computational forward model is replaced by a large dictionary of possible outcomes, and the problem is to identify the dictionary entries that best match the data, akin to traditional query matching in search engines. Sparse coding techniques are used to guarantee that the dictionary matching identifies only few of the dictionary entries, and dictionary compression methods are used to reduce the complexity of the matching problem. In this article, we propose a work flow to facilitate the dictionary matching process. First, the full dictionary is divided into subdictionaries that are separately compressed. The error introduced by the dictionary compression is handled in the Bayesian framework as a modeling error. Furthermore, we propose a new Bayesian data-driven group sparsity coding method to help identify subdictionaries that are not relevant for the dictionary matching. After discarding irrelevant subdictionaries, the dictionary matching is addressed as a deflated problem using sparse coding. The compression and deflation steps can lead to substantial decreases of the computational complexity. The effectiveness of compensating for the dictionary compression error and using the novel group sparsity promotion to deflate the original dictionary are illustrated by applying the methodology to real world problems, the glitch detection in the LIGO experiment and hyperspectral remote sensing.

Bayesian hierarchical models have been demonstrated to provide efficient algorithms for finding sparse solutions to ill-posed inverse problems. The models comprise typically a conditionally Gaussian prior model for the unknown, augmented by a hyperprior model for the variances. A widely used choice for the hyperprior is a member of the family of generalized gamma distributions. Most of the work in the literature has concentrated on numerical approximation of the maximum a posteriori (MAP) estimates, and less attention has been paid on sampling methods or other means for uncertainty quantification. Sampling from the hierarchical models is challenging mainly for two reasons: The hierarchical models are typically high-dimensional, thus suffering from the curse of dimensionality, and the strong correlation between the unknown of interest and its variance can make sampling rather inefficient. This work addresses mainly the first one of these obstacles. By using a novel reparametrization, it is shown how the posterior distribution can be transformed into one dominated by a Gaussian white noise, allowing sampling by using the preconditioned Crank-Nicholson (pCN) scheme that has been shown to be efficient for sampling from distributions dominated by a Gaussian component. Furthermore, a novel idea for speeding up the pCN in a special case is developed, and the question of how strongly the hierarchical models are concentrated on sparse solutions is addressed in light of a computed example.

Stellar streams are sensitive probes of the Galactic potential. The likelihood of a stream model given stream data is often assessed using simulations. However, comparing to simulations is challenging when even the stream paths can be hard to quantify. Here we present a novel application of Self-Organizing Maps and first-order Kalman Filters to reconstruct a stream's path, propagating measurement errors and data sparsity into the stream path uncertainty. The technique is Galactic-model independent, non-parametric, and works on phase-wrapped streams. With this technique, we can uniformly analyze and compare data with simulations, enabling both comparison of simulation techniques and ensemble analysis with stream tracks of many stellar streams. Our method is implemented in the public Python package TrackStream, available at https://github.com/nstarman/trackstream., Comment: 16 pages, 11 figures, preprint

The different active roles of neurons and astrocytes during neuronal activation are associated with the metabolic processes necessary to supply the energy needed for their respective tasks at rest and during neuronal activation. Metabolism, in turn, relies on the delivery of metabolites and removal of toxic byproducts through diffusion processes and the cerebral blood flow. A comprehensive mathematical model of brain metabolism should account not only for the biochemical processes and the interaction of neurons and astrocytes, but also the diffusion of metabolites. In the present article, we present a computational methodology based on a multidomain model of the brain tissue and a homogenization argument for the diffusion processes. In our spatially distributed compartment model, communication between compartments occur both through local transport fluxes, as is the case within local astrocyte-neuron complexes, and through diffusion of some substances in some of the compartments. The model assumes that diffusion takes place in the extracellular space (ECS) and in the astrocyte compartment. In the astrocyte compartment, the diffusion across the syncytium network is implemented as a function of gap junction strength. The diffusion process is implemented numerically by means of a finite element method (FEM) based spatial discretization, and robust stiff solvers are used to time integrate the resulting large system. Computed experiments show the effects of ECS tortuosity, gap junction strength and spatial anisotropy in the astrocyte network on the brain energy metabolism.

A common task in inverse problems and imaging is finding a solution that is sparse, in the sense that most of its components vanish. In the framework of compressed sensing, general results guaranteeing exact recovery have been proven. In practice, sparse solutions are often computed combining $\ell_1$-penalized least squares optimization with an appropriate numerical scheme to accomplish the task. A computationally efficient alternative for finding sparse solutions to linear inverse problems is provided by Bayesian hierarchical models, in which the sparsity is encoded by defining a conditionally Gaussian prior model with the prior parameter obeying a generalized gamma distribution. An iterative alternating sequential (IAS) algorithm has been demonstrated to lead to a computationally efficient scheme, and combined with Krylov subspace iterations with an early termination condition, the approach is particularly well suited for large scale problems. Here the Bayesian approach to sparsity is extended to problems whose solution allows a sparse coding in an overcomplete system such as composite frames. It is shown that among the multiple possible representations of the unknown, the IAS algorithm, and in particular, a hybrid version of it, is effectively identifying the most sparse solution. Computed examples show that the method is particularly well suited not only for traditional imaging applications but also for dictionary learning problems in the framework of machine learning.

Mathematical models of SARS-CoV-2 spread are used for guiding the design of mitigation steps aimed at containing and decelerating the contagion, and at identifying impending breaches of health care system surge capacity. The challenges of having only lacunary information about daily new infections are compounded by the geographic heterogeneity of the population. To address this problem, we propose to account for the differences between rural and urban settings using network-based, distributed models where the spread of the pandemic is described in distinct local cohorts with nested SEIR models. The setting of the model parameters takes into account the fact that SARS-CoV-2 transmission occurs mostly via human-to-human contact, and that the frequency of contact among individuals differs between urban and rural areas, and may change over time. Moreover, the probability that the virus spreads into an uninfected community is associated with influx of individuals from other communities where the infection is present. To account for these important aspects, each node of the network is characterized by the frequency of contact between its members and by its level of connectivity with other nodes. Census and cell phone data can be used to set up the adjacency matrix of the network, which can, in turn, be modified to account for different levels of mitigation measures. In order to make the network SEIR model that we propose easy to customize, it is formulated in terms of easily interpretable parameters that can be estimated from available community level data. The models parameters are estimated with Bayesian techniques using COVID-19 data for the states of Ohio and Michigan. The network model also gives rise to a geographically distributed computational model that explains the geographic dynamics of the contagion, e.g., in larger cities surrounded by suburban and rural areas., Comment: 21 pages, 15 figures

In this article, we consider a dynamic epidemiology model for the spread of the COVID-19 infection. Starting from the classical SEIR model, the model is modified so as to better describe characteristic features of the underlying pathogen and its infectious modes. In line with the large number of secondary infections not related to contact with documented infectious individuals, the model includes a cohort of asymptomatic or oligosymptomatic infectious individuals, not accounted for in the data of new daily counts of infections. A Bayesian particle filtering algorithm is used to update dynamically the relevant cohort and simultaneously estimate the transmission rate as the new data on the number of new infections and disease related death become available. The underlying assumption of the model is that the infectivity rate is dynamically changing during the epidemics, either because of a mutation of the pathogen or in response to mitigation and containment measures. The sequential Bayesian framework naturally provides a quantification of the uncertainty in the estimate of the model parameters, including the reproduction number, and of the size of the different cohorts. Moreover, we introduce a dimensionless quantity, which is the equilibrium ratio between asymptomatic and symptomatic cohort sizes, and propose a simple formula to estimate the quantity. This ratio leads naturally to another dimensionless quantity that plays the role of the basic reproduction number $R_0$ of the model. When we apply the model and particle filter algorithm to COVID-19 infection data from several counties in Northeastern Ohio and Southeastern Michigan we found the proposed reproduction number $R_0$ to have a consistent dynamic behavior within both states, thus proving to be a reliable summary of the success of the mitigation measures., Comment: 21 pages, 8 figures

The recovery of sparse generative models from few noisy measurements is an important and challenging problem. Many deterministic algorithms rely on some form of $\ell_1$-$\ell_2$ minimization to combine the computational convenience of the $\ell_2$ penalty and the sparsity promotion of the $\ell_1$. It was recently shown within the Bayesian framework that sparsity promotion and computational efficiency can be attained with hierarchical models with conditionally Gaussian priors and gamma hyperpriors. The related Gibbs energy function is a convex functional and its minimizer, which is the MAP estimate of the posterior, can be computed efficiently with the globally convergent Iterated Alternating Sequential (IAS) algorithm \cite{CSS}. Generalization of the hyperpriors for these sparsity promoting hierarchical models to generalized gamma family yield either globally convex Gibbs energy functionals, or can exhibit local convexity for some choices for the hyperparameters. \cite{CPrSS}. The main problem in computing the MAP solution for greedy hyperpriors that strongly promote sparsity is the presence of local minima. To overcome the premature stopping at a spurious local minimizer, we propose two hybrid algorithms that first exploit the global convergence associated with gamma hyperpriors to arrive in a neighborhood of the unique minimizer, then adopt a generalized gamma hyperprior that promote sparsity more strongly. The performance of the two algorithms is illustrated with computed examples.

A recently proposed IAS MEG inverse solver algorithm, based on the coupling of a hierarchical Bayesian model with computationally efficient Krylov subspace linear solver, has been shown to perform well for both superficial and deep brain sources. However, a systematic study of its sensitivity and specificity as a function of the activity location is still missing. We propose novel statistical protocols to quantify the performance of MEG inverse solvers, focusing in particular on their sensitivity and specificity in identifying active brain regions. We use these protocols for a systematic study of the sensitivity and specificity of the IAS MEG inverse solver, comparing the performance with three standard inversion methods, wMNE, dSPM, and sLORETA. To avoid the bias of anecdotal tests towards a particular algorithm, the proposed protocols are Monte Carlo sampling based, generating an ensemble of activity patches in each brain region identified in a given atlas. The sensitivity is measured by how much, on average, the reconstructed activity is concentrated in the brain region of the simulated active patch. The specificity analysis is based on Bayes factors, interpreting the estimated current activity as data for testing the hypothesis that the active brain region is correctly identified, vs. the hypothesis of any erroneous attribution. The methodology allows the presence of a single or several simultaneous activity regions, without assuming the knowledge of the number of active regions. The testing protocols suggest that the IAS solver performs well in terms of sensitivity and specificity both with cortical and subcortical activity estimation.

In computational inverse problems, it is common that a detailed and accurate forward model is approximated by a computationally less challenging substitute. The model reduction may be necessary to meet constraints in computing time when optimization algorithms are used to find a single estimate, or to speed up Markov chain Monte Carlo (MCMC) calculations in the Bayesian framework. The use of an approximate model introduces a discrepancy, or modeling error, that may have a detrimental effect on the solution of the ill-posed inverse problem, or it may severely distort the estimate of the posterior distribution. In the Bayesian paradigm, the modeling error can be considered as a random variable, and by using an estimate of the probability distribution of the unknown, one may estimate the probability distribution of the modeling error and incorporate it into the inversion. We introduce an algorithm which iterates this idea to update the distribution of the model error, leading to a sequence of posterior distributions that are demonstrated empirically to capture the underlying truth with increasing accuracy. Since the algorithm is not based on rejections, it requires only limited full model evaluations. We show analytically that, in the linear Gaussian case, the algorithm converges geometrically fast with respect to the number of iterations. For more general models, we introduce particle approximations of the iteratively generated sequence of distributions; we also prove that each element of the sequence converges in the large particle limit. We show numerically that, as in the linear case, rapid convergence occurs with respect to the number of iterations. Additionally, we show through computed examples that point estimates obtained from this iterative algorithm are superior to those obtained by neglecting the model error., Comment: 39 pages, 9 figures

The Linear Multistep Method Particle Filter (LMM PF) is a method for predicting the evolution in time of a evolutionary system governed by a system of differential equations. If some of the parameters of the governing equations are unknowns, it is possible to organize the calculations so as to estimate them while following the evolution of the system in time. The underlying assumption in the approach that we present is that all unknowns are modelled as random variables, where the randomness is an indication of the uncertainty of their values rather than an intrinsic property of the quantities. Consequently, the states of the system and the parameters are described in probabilistic terms by their density, often in the form of representative samples. This approach is particularly attractive in the context of parameter estimation inverse problems, because the statistical formulation naturally provides a means of assessing the uncertainty in the solution via the spread of the distribution. The computational efficiency of the underlying sampling technique is crucial for the success of the method, because the accuracy of the solution depends on the ability to produce representative samples from the distribution of the unknown parameters. In this paper LMM PF is tested on a skeletal muscle metabolism problem, which was previously treated within the Ensemble Kalman filtering framework. Here numerical evidences are used to highlight the correlation between the main sources of errors and the influence of the linera multistep method adopted. Finally, we analyzed the effect of replacing LMM with Runge-Kutta class integration methods for supporting the PF technique.