Erkki Somersalo

A way to lower computational cost in large scale inverse problems and problems depending on poorly known model parameters is to replace the detailed model by an approximate one. Inverse problems are typically ill-posed, and the model discrepancy introduced by using approximate models often shows up in the computed solutions as disturbing artifacts or blurring. In this article, we consider two methods of addressing certain types of modeling errors, the Bayesian approximation error (BAE) method and linear algebraic spotlight inversion to suppress clutter in the computational model by orthogonal projections. Through the process of analyzing the two approaches, we show that they turn out to be closely related but not equivalent, and we highlight a connection to sketching schemes in randomized linear algebra. The similarities between the methods and their successful suppression of most of the clutter effects is elucidated with two computed examples, one addressing of x-ray tomography and the other electrical impedance tomography.

Many computational problems involve solving a linear system of equations, although only a subset of the entries of the solution are needed. In inverse problems, where the goal is to estimate unknown parameters from indirect noisy observations, it is not uncommon that the forward model linking the observed variables to the unknowns depends on variables that are not of primary interest, often referred to as nuisance parameters. In this article, we consider linear problems and propose a novel projection technique to eliminate, or at least mitigate, the contribution of the nuisance parameters in the model. We refer to this approach as spotlight inversion, as it allows to focus on only the portion of primary interest of the unknown parameter vector, leaving the uninteresting part in the shadow. The viability of the approach is illustrated with two computed examples, one where it works as model reduction for a finite element approximation of an elliptic PDE, the other amounting to local fanbeam X-ray tomography, spotlighting the region of interest that is part of the full target.

A central goal in many brain studies is the identification of those brain regions that are activated during an observation window that may correspond to a motor task, a stimulus, or simply a resting state. While functional MRI is currently the most commonly employed modality for such task, methods based on the electromagnetic activity of the brain are valuable alternatives because of their excellent time resolution and of the fact that the measured signals are directly related to brain activation and not to a secondary effect such as the hemodynamic response. In this work we focus on the MEG modality, investigating the performance of a recently proposed Bayesian dictionary learning (BDL) algorithm for brain region identification. The partitioning of the source space into the 148 regions of interest (ROI) corresponding to parcellation of the Destrieux atlas provides a natural determination of the subdictionaries necessary for the BDL algorithm. We design a simulation protocol where a small randomly selected patch in each ROI is activated, the MEG signal is computed and the inverse problem of active brain region identification is solved using the BDL algorithm. The BDL algorithm consists of two phases, the first one comprising dictionary compression and Bayesian compression error analysis, and the second one performing dictionary coding with a deflated dictionary built on the output of the first phase, both steps relying on Bayesian sparsity promoting computations. For assessing the performance, we give a probabilistic interpretation of the confusion matrix, and consider different impurity measures for a multi-class classifier. Key words: Sparse coding, Bayesian approximation error, Matrix factorization, Impurity index, Bayesian hierarchical models

Gas transport across cell membrane is a very important process in cell biology, as it is essential for many crucial tasks, including cell respiration and pH regulation. Since it is difficult to observe directly how gases cross cell membranes, mathematical models have played and continue to play an important role in providing a context for the experimental findings depending on gas exchange, and to test different hypotheses about the gas passage mechanism. The first family of gas transport models has been able to match qualitatively the experimental data, and a second generation of computational models has reproduced and explained the dynamic range of the data. These models depend on parameters that cannot be measured experimentally, but can be estimated from the available data and some a priori knowledge about the problem based on first principles of cell physiology. The membrane permeability is arguably one of the most important parameters and can also be used also to shed light on the validity of different models of the transport mechanisms. In this contribution we propose a new, computationally efficient approach to the cell permeability estimation from measurements of the pH on the cell membrane over the course of a lab experiment. The methodology is inspired by a technique used in data science applications known as dictionary learning. Relevance to Life Sciences.: In the late 1990s, the suggestion that gases are transported via preferred gas channels embedded into the cell membrane challenged the century-old Overton theory that gases pass through the lipid cell membrane by diffusing across it driven by the concentration gradient. Since the laboratory experiments alone do not provide enough evidence to favor one or another of the proposed mechanisms, mathematical models have been introduced to provide a context for the interpretation of the measurements. In this article we address the problem of estimating intrinsic properties of the cell membrane that could help resolve the question about the role of aquaporins and other membrane-embedded proteins in gas transport across a cell membrane from typical pH measurements on the membrane of the oocytes of Xenopus laevis. Mathematical Content.: Following up on a previous work where the membrane permeability was estimated using particle filter techniques, in this article we propose an algorithm based on dictionary learning for estimating cell membrane permeability and other pertinent model parameters. Dictionary learning methods have proven to provide versatile tools to solve inverse problems in which the insufficiency of the available data or the complexity of the forward model makes the use of optimization-based inverse solvers impractical. Computed examples illustrate that the novel approach, which can be applied when the properties of the membrane do not change in the course of the data collection process, is computationally more efficient than the particle filter approach.

It is well-known that the posterior density of linear inverse problems with Gaussian prior and Gaussian likelihood is also Gaussian, hence completely described by its covariance and expectation. Sampling from a Gaussian posterior may be important in the analysis of various non-Gaussian inverse problems in which estimates from a Gaussian posterior distribution constitute an intermediate stage in a Bayesian workflow. Sampling from a Gaussian distribution is straightforward if the Cholesky factorization of the covariance matrix or its inverse is available; however, when the unknown is high dimensional, the computation of the posterior covariance may be unfeasible. If the linear inverse problem is underdetermined, it is possible to exploit the orthogonality of the fundamental subspaces associated with the coefficient matrix together with the idea behind the randomize-then-optimize approach to design a low complexity posterior sampler that does not require the posterior covariance to be formed. The performance of the proposed sampler is illustrated with a few computed examples, including non-Gaussian problems with nonlinear forward model, and hierarchical models comprising a conditionally Gaussian submodel.

Optical flow velocimetry (OFV) is a method for determining dense and accurate velocity fields from a pair of particle images by solving the classical optical flow problem. However, this is an ill-posed inverse problem, which generally entails minimizing a weighted sum of two terms–fidelity and regularization–and the weights in the sum are parameters that require manual tuning based on the properties of both the flow and the particle images. This manual tuning has historically been a consistent challenge that has limited the general applicability of OFV for experimental data, as the calculated velocity field is sensitive to the value of the weights. This work proposes a hierarchical model for the weighting parameters in the framework of a maximum a posteriori-based Bayesian optimization approach. The method replaces the classical Lagrange multiplier weighting parameter with a new, less-sensitive parameter that can be automatically predetermined from experimental images. The resulting method is tested on three different synthetic particle image velocimetry (PIV) datasets and on experimental particle images. The method is found to be capable of self-adjusting the local weights of the optimization process in real-time while simultaneously determining the velocity field, leading to an optimally regularized estimate of the velocity field without requiring any dataset-specific manual tuning of the parameters. The presented approach is the first truly general, parameter-free optical flow method for particle image velocimetry (PIV) images. The developed method is freely available as a part of the PIVlab package.

The inherent redundancy of the musculoskeletal systems is traditionally solved by optimizing a cost function. This approach may not be correct to model non-adult or pathological populations likely to adopt a “non-optimal” motor control strategy. Over the years, various methods have been developed to address this limitation, such as the stochastic approach. A well-known implementation of this approach, Metabolica, samples a wide number of plausible solutions instead of searching for a single one, leveraging Bayesian statistics and Markov Chain Monte Carlo algorithm, yet allowing muscles to abruptly change their activation levels. To overcome this and other limitations, we developed a new implementation of the stochastic approach (Myobolica), adding constraints and parameters to ensure the identification of physiological solutions. The aim of this study was to evaluate Myobolica, and quantify the differences in terms of width of the solution band (muscle control variability) compared to Metabolica. To this end, both muscle forces and knee joint force solutions bands estimated by the two approaches were compared to one another, and against (i) the solution identified by static optimization and (ii) experimentally measured knee joint forces. The use of Myobolica led to a marked narrowing of the solution band compared to Metabolica. Furthermore, the Myobolica solutions well correlated with the experimental data (R ² = 0.92, RMSE = 0.3 BW), but not as much with the optimal solution (R ² = 0.82, RMSE = 0.63 BW). Additional analyses are required to confirm the findings and further improve this implementation.

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 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.

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.