Bangti Jin

Bangti Jin, Raytcho Lazarov, Zhi Zhou

The subdiffusion equation with a Caputo fractional derivative of order $α\in(0,1)$ in time arises in a wide variety of practical applications, and it is often adopted to model anomalous subdiffusion processes in heterogeneous media. The L1 scheme is one of the most popular and successful numerical methods for discretizing the Caputo fractional derivative in time. The scheme was analyzed earlier independently by Lin and Xu (2007) and Sun and Wu (2006), and an $O(τ^{2-α})$ convergence rate was established, under the assumption that the solution is twice continuously differentiable in time. However, in view of the smoothing property of the subdiffusion equation, this regularity condition is restrictive, since it does not hold even for the homogeneous problem with a smooth initial data. In this work, we revisit the error analysis of the scheme, and establish an $O(τ)$ convergence rate for both smooth and nonsmooth initial data. The analysis is valid for more general sectorial operators. In particular, the L1 scheme is applied to one-dimensional space-time fractional diffusion equations, which involves also a Riemann-Liouville derivative of order $β\in(3/2,2)$ in space, and error estimates are provided for the fully discrete scheme. Numerical experiments are provided to verify the sharpness of the error estimates, and robustness of the scheme with respect to data regularity.

Bangti Jin, Raytcho Lazarov, Zhi Zhou

We consider the initial boundary value problem for the homogeneous time-fractional diffusion equation $\partial^α_t u - \De u =0$ ($0< α< 1$) with initial condition $u(x,0)=v(x)$ and a homogeneous Dirichlet boundary condition in a bounded polygonal domain $Ω$. We shall study two semidiscrete approximation schemes, i.e., Galerkin FEM and lumped mass Galerkin FEM, by using piecewise linear functions. We establish optimal with respect to the regularity of the solution error estimates, including the case of nonsmooth initial data, i.e., $v \in L_2(Ω)$.

Bangti Jin, Buyang Li, Zhi Zhou

We develop proper correction formulas at the starting $k-1$ steps to restore the desired $k^{\rm th}$-order convergence rate of the $k$-step BDF convolution quadrature for discretizing evolution equations involving a fractional-order derivative in time. The desired $k^{\rm th}$-order convergence rate can be achieved even if the source term is not compatible with the initial data, which is allowed to be nonsmooth. We provide complete error estimates for the subdiffusion case $α\in (0,1)$, and sketch the proof for the diffusion-wave case $α\in(1,2)$. Extensive numerical examples are provided to illustrate the effectiveness of the proposed scheme.

Bangti Jin, Buyang Li, Zhi Zhou

We present a general framework for the rigorous numerical analysis of time-fractional nonlinear parabolic partial differential equations, with a fractional derivative of order $α\in(0,1)$ in time. The framework relies on three technical tools: a fractional version of the discrete Grönwall-type inequality, discrete maximal regularity, and regularity theory of nonlinear equations. We establish a general criterion for showing the fractional discrete Grönwall inequality, and verify it for the L1 scheme and convolution quadrature generated by BDFs. Further, we provide a complete solution theory, e.g., existence, uniqueness and regularity, for a time-fractional diffusion equation with a Lipschitz nonlinear source term. Together with the known results of discrete maximal regularity, we derive pointwise $L^2(Ω)$ norm error estimates for semidiscrete Galerkin finite element solutions and fully discrete solutions, which are of order $O(h^2)$ (up to a logarithmic factor) and $O(τ^α)$, respectively, without any extra regularity assumption on the solution or compatibility condition on the problem data. The sharpness of the convergence rates is supported by the numerical experiments.

Emilia Bazhlekova, Bangti Jin, Raytcho Lazarov et al.

We study the Rayleigh-Stokes problem for a generalized second-grade fluid which involves a Riemann-Liouville fractional derivative in time, and present an analysis of the problem in the continuous, space semidiscrete and fully discrete formulations. We establish the Sobolev regularity of the homogeneous problem for both smooth and nonsmooth initial data $v$, including $v\in L^2(Ω)$. A space semidiscrete Galerkin scheme using continuous piecewise linear finite elements is developed, and optimal with respect to initial data regularity error estimates for the finite element approximations are derived. Further, two fully discrete schemes based on the backward Euler method and second-order backward difference method and the related convolution quadrature are developed, and optimal error estimates are derived for the fully discrete approximations for both smooth and nonsmooth initial data. Numerical results for one- and two-dimensional examples with smooth and nonsmooth initial data are presented to illustrate the efficiency of the method, and to verify the convergence theory.

Imraj RD Singh, Alexander Denker, Riccardo Barbano et al.

Score-based generative models have demonstrated highly promising results for medical image reconstruction tasks in magnetic resonance imaging or computed tomography. However, their application to Positron Emission Tomography (PET) is still largely unexplored. PET image reconstruction involves a variety of challenges, including Poisson noise with high variance and a wide dynamic range. To address these challenges, we propose several PET-specific adaptations of score-based generative models. The proposed framework is developed for both 2D and 3D PET. In addition, we provide an extension to guided reconstruction using magnetic resonance images. We validate the approach through extensive 2D and 3D $\textit{in-silico}$ experiments with a model trained on patient-realistic data without lesions, and evaluate on data without lesions as well as out-of-distribution data with lesions. This demonstrates the proposed method's robustness and significant potential for improved PET reconstruction.

Riccardo Barbano, Javier Antorán, Johannes Leuschner et al. · cambridge

Deep learning has been widely used for solving image reconstruction tasks but its deployability has been held back due to the shortage of high-quality training data. Unsupervised learning methods, such as the deep image prior (DIP), naturally fill this gap, but bring a host of new issues: the susceptibility to overfitting due to a lack of robust early stopping strategies and unstable convergence. We present a novel approach to tackle these issues by restricting DIP optimisation to a sparse linear subspace of its parameters, employing a synergy of dimensionality reduction techniques and second order optimisation methods. The low-dimensionality of the subspace reduces DIP's tendency to fit noise and allows the use of stable second order optimisation methods, e.g., natural gradient descent or L-BFGS. Experiments across both image restoration and tomographic tasks of different geometry and ill-posedness show that second order optimisation within a low-dimensional subspace is favourable in terms of optimisation stability to reconstruction fidelity trade-off.

Bangti Jin, Buyang Li, Zhi Zhou

In this work, we analyze a Crank-Nicolson type time stepping scheme for the subdiffusion equation, which involves a Caputo fractional derivative of order $α\in (0,1)$ in time. It hybridizes the backward Euler convolution quadrature with a $θ$-type method, with the parameter $θ$ dependent on the fractional order $α$ by $θ=α/2$, and naturally generalizes the classical Crank-Nicolson method. We develop essential initial corrections at the starting two steps for the Crank-Nicolson scheme, and together with the Galerkin finite element method in space, obtain a fully discrete scheme. The overall scheme is easy to implement, and robust with respect to data regularity. A complete error analysis of the fully discrete scheme is provided, and a second-order accuracy in time is established for both smooth and nonsmooth problem data. Extensive numerical experiments are provided to illustrate its accuracy, efficiency and robustness, and a comparative study also indicates its competitive with existing schemes.

Bangti Jin, Buyang Li, Zhi Zhou

In this work, we establish the maximal $\ell^p$-regularity for several time stepping schemes for a fractional evolution model, which involves a fractional derivative of order $α\in(0,2)$, $α\neq 1$, in time. These schemes include convolution quadratures generated by backward Euler method and second-order backward difference formula, the L1 scheme, explicit Euler method and a fractional variant of the Crank-Nicolson method. The main tools for the analysis include operator-valued Fourier multiplier theorem due to Weis [48] and its discrete analogue due to Blunck [10]. These results generalize the corresponding results for parabolic problems.

Bangti Jin, Raytcho Lazarov, Zhi Zhou

In this work, we develop variational formulations of Petrov-Galerkin type for one-dimensional fractional boundary value problems involving either a Riemann-Liouville or Caputo derivative of order $α\in(3/2, 2)$ in the leading term and both convection and potential terms. They arise in the mathematical modeling of asymmetric super-diffusion processes in heterogeneous media. The well-posedness of the formulations and sharp regularity pickup of the variational solutions are established. A novel finite element method is developed, which employs continuous piecewise linear finite elements and "shifted" fractional powers for the trial and test space, respectively. The new approach has a number of distinct features: It allows deriving optimal error estimates in both $L^2(D)$ and $H^1(D)$ norms; and on a uniform mesh, the stiffness matrix of the leading term is diagonal and the resulting linear system is well conditioned. Further, in the Riemann-Liouville case, an enriched FEM is proposed to improve the convergence. Extensive numerical results are presented to verify the theoretical analysis and robustness of the numerical scheme.

Kazufumi Ito, Bangti Jin, Tomoya Takeuchi

We study multi-parameter Tikhonov regularization, i.e., with multiple penalties. Such models are useful when the sought-for solution exhibits several distinct features simultaneously. Two choice rules, i.e., discrepancy principle and balancing principle, are studied for choosing an appropriate (vector-valued) regularization parameter, and some theoretical results are presented. In particular, the consistency of the discrepancy principle as well as convergence rate are established, and an a posteriori error estimate for the balancing principle is established. Also two fixed point algorithms are proposed for computing the regularization parameter by the latter rule. Numerical results for several nonsmooth multi-parameter models are presented, which show clearly their superior performance over their single-parameter counterparts.

Bangti Jin, Buyang Li, Zhi Zhou

In this work, we present numerical analysis for a distributed optimal control problem, with box constraint on the control, governed by a subdiffusion equation which involves a fractional derivative of order $α\in(0,1)$ in time. The fully discrete scheme is obtained by applying the conforming linear Galerkin finite element method in space, L1 scheme/backward Euler convolution quadrature in time, and the control variable by a variational type discretization. With a space mesh size $h$ and time stepsize $τ$, we establish the following order of convergence for the numerical solutions of the optimal control problem: $O(τ^{\min({1}/{2}+α-ε,1)}+h^2)$ in the discrete $L^2(0,T;L^2(Ω))$ norm and $O(τ^{α-ε}+\ell_h^2h^2)$ in the discrete $L^\infty(0,T;L^2(Ω))$ norm, with any small $ε>0$ and $\ell_h=\ln(2+1/h)$. The analysis relies essentially on the maximal $L^p$-regularity and its discrete analogue for the subdiffusion problem. Numerical experiments are provided to support the theoretical results.

Tobias Kluth, Bangti Jin

Magnetic particle imaging (MPI) is a medical imaging modality of recent origin, and it exploits the nonlinear magnetization phenomenon to recover the spatially dependent concentration of the nanoparticles. Currently, image reconstruction in MPI is frequently carried out by standard Tikhonov regularization with nonnegativity constraint, which is then minimized by a Kaczmarz type method. In this work, we revisit several issues in the numerical reconstruction in MPI from the perspective of modern inverse theory, i.e., the choice of data fidelity, and choosing a suitable regularization parameter and accelerating Kaczmarz iteration via randomized singular value decomposition. These algorithmic tricks are straightforward to implement and easy to incorporate in existing reconstruction algorithms. Their significant potentials are illustrated by extensive numerical experiments on a publicly available dataset.

Bangti Jin, Raytcho Lazarov, Xiliang Lu et al.

We consider a boundary value problem involving a Riemann-Liouville fractional derivative of order $α\in (3/2,2)$ on the unit interval $(0,1)$. The standard Galerkin finite element approximation converges slowly due to the presence of singularity term $x^{α-1}$ in the solution representation. In this work, we develop a simple technique, by transforming it into a second-order two-point boundary value problem with nonlocal low order terms, whose solution can reconstruct directly the solution to the original problem. The stability of the variational formulation, and the optimal regularity pickup of the solution are analyzed. A novel Galerkin finite element method with piecewise linear or quadratic finite elements is developed, and $L^2(D)$ error estimates are provided. The approach is then applied to the corresponding fractional Sturm-Liouville problem, and error estimates of the eigenvalue approximations are given. Extensive numerical results fully confirm our theoretical study.

Bangti Jin, Zhi Zhou

In this work, we develop a novel Galerkin-L1-POD scheme for the subdiffusion model with a Caputo fractional derivative of order $α\in (0,1)$ in time, which is often used to describe anomalous diffusion processes in heterogeneous media. The nonlocality of the fractional derivative requires storing all the solutions from time zero. The proposed scheme is based on continuous piecewise linear finite elements, L1 time stepping, and proper orthogonal decomposition (POD). By constructing an effective reduced-order scheme using problem-adapted basis functions, it can significantly reduce the computational complexity and storage requirement. We shall provide a complete error analysis of the scheme under realistic regularity assumptions by means of a novel energy argument. Extensive numerical experiments are presented to verify the convergence analysis and the efficiency of the proposed scheme.

Bangti Jin, Zhi Zhou

We consider a two-point boundary value problem involving a Riemann-Liouville fractional derivative of order $\al\in (1,2)$ in the leading term on the unit interval $(0,1)$. Generally the standard Galerkin finite element method can only give a low-order convergence even if the source term is very smooth due to the presence of the singularity term $x^{\al-1}$ in the solution representation. In order to enhance the convergence, we develop a simple singularity reconstruction strategy by splitting the solution into a singular part and a regular part, where the former captures explicitly the singularity. We derive a new variational formulation for the regular part, and establish that the Galerkin approximation of the regular part can achieve a better convergence order in the $L^2(0,1)$, $H^{\al/2}(0,1)$ and $L^\infty(0,1)$-norms than the standard Galerkin approach, with a convergence rate for the recovered singularity strength identical with the $L^2(0,1)$ error estimate. The reconstruction approach is very flexible in handling explicit singularity, and it is further extended to the case of a Neumann type boundary condition on the left end point, which involves a strong singularity $x^{\al-2}$. Extensive numerical results confirm the theoretical study and efficiency of the proposed approach.

Siyu Cen, Bangti Jin, Kwancheol Shin et al.

Electrical impedance tomography (EIT) is a noninvasive medical imaging modality utilizing the current-density/voltage data measured on the surface of the subject. Calderón's method is a relatively recent EIT imaging algorithm that is non-iterative, fast, and capable of reconstructing complex-valued electric impedances. However, due to the regularization via low-pass filtering and linearization, the reconstructed images suffer from severe blurring and under-estimation of the exact conductivity values. In this work, we develop an enhanced version of Calderón's method, using {deep} convolution neural networks (i.e., U-net) {as an effective targeted post-processing step, and term the resulting method by deep Calderón's method.} Specifically, we learn a U-net to postprocess the EIT images generated by Calderón's method so as to have better resolutions and more accurate estimates of conductivity values. We simulate chest configurations with which we generate the current-density/voltage boundary measurements and the corresponding reconstructed images by Calderón's method. With the paired training data, we learn the deep neural network and evaluate its performance on real tank measurement data. The experimental results indicate that the proposed approach indeed provides a fast and direct (complex-valued) impedance tomography imaging technique, and substantially improves the capability of the standard Calderón's method.

Bangti Jin, Raytcho Lazarov, Joseph Pasciak et al.

We investigate semi-discrete numerical schemes based on the standard Galerkin and lumped mass Galerkin finite element methods for an initial-boundary value problem for homogeneous fractional diffusion problems with non-smooth initial data. We assume that $Ω\subset \mathbb{R}^d$, $d=1,2,3$ is a convex polygonal (polyhedral) domain. We theoretically justify optimal order error estimates in $L_2$- and $H^1$-norms for initial data in $H^{-s}(Ω),~0\le s \le 1$. We confirm our theoretical findings with a number of numerical tests that include initial data $v$ being a Dirac $δ$-function supported on a $(d-1)$-dimensional manifold.

Bangti Jin, Xiyao Li, Xiliang Lu

Conductivity imaging represents one of the most important tasks in medical imaging. In this work we develop a neural network based reconstruction technique for imaging the conductivity from the magnitude of the internal current density. It is achieved by formulating the problem as a relaxed weighted least-gradient problem, and then approximating its minimizer by standard fully connected feedforward neural networks. We derive bounds on two components of the generalization error, i.e., approximation error and statistical error, explicitly in terms of properties of the neural networks (e.g., depth, total number of parameters, and the bound of the network parameters). We illustrate the performance and distinct features of the approach on several numerical experiments. Numerically, it is observed that the approach enjoys remarkable robustness with respect to the presence of data noise.

Riccardo Barbano, Johannes Leuschner, Javier Antorán et al.

We investigate adaptive design based on a single sparse pilot scan for generating effective scanning strategies for computed tomography reconstruction. We propose a novel approach using the linearised deep image prior. It allows incorporating information from the pilot measurements into the angle selection criteria, while maintaining the tractability of a conjugate Gaussian-linear model. On a synthetically generated dataset with preferential directions, linearised DIP design allows reducing the number of scans by up to 30% relative to an equidistant angle baseline.

Victor M. Calo, Nathan Collier, Matthias Gehre et al.

We study the numerical recovery of Manning's roughness coefficient for the diffusive wave approximation of the shallow water equation. We describe a conjugate gradient method for the numerical inversion. Numerical results for one-dimensional model are presented to illustrate the feasibility of the approach. Also we provide a proof of the differentiability of the weak form with respect to the coefficient as well as the continuity and boundedness of the linearized operator under reasonable assumptions using the maximal parabolic regularity theory.

Yongcheng Dai, Bangti Jin, Ramesh Sau et al.

In this work, we investigate a neural network based solver for optimal control problems (without / with box constraint) for linear and semilinear second-order elliptic problems. It utilizes a coupled system derived from the first-order optimality system of the optimal control problem, and employs deep neural networks to represent the solutions to the reduced system. We present an error analysis of the scheme, and provide $L^2(Ω)$ error bounds on the state, control and adjoint in terms of neural network parameters (e.g., depth, width, and parameter bounds) and the numbers of sampling points. The main tools in the analysis include offset Rademacher complexity and boundedness and Lipschitz continuity of neural network functions. We present several numerical examples to illustrate the method and compare it with two existing ones.

Tianhao Hu, Bangti Jin, Zhi Zhou

In this work, we develop an efficient solver based on neural networks for second-order elliptic equations with variable coefficients and singular sources. This class of problems covers general point sources, line sources and the combination of point-line sources, and has a broad range of practical applications. The proposed approach is based on decomposing the true solution into a singular part that is known analytically using the fundamental solution of the Laplace equation and a regular part that satisfies a suitable modified elliptic PDE with a smoother source, and then solving for the regular part using the deep Ritz method. A path-following strategy is suggested to select the penalty parameter for enforcing the Dirichlet boundary condition. Extensive numerical experiments in two- and multi-dimensional spaces with point sources, line sources or their combinations are presented to illustrate the efficiency of the proposed approach, and a comparative study with several existing approaches based on neural networks is also given, which shows clearly its competitiveness for the specific class of problems. In addition, we briefly discuss the error analysis of the approach.

Bangti Jin, Xiyao Li, Qimeng Quan et al.

In this work we develop a novel approach using deep neural networks to reconstruct the conductivity distribution in elliptic problems from one measurement of the solution over the whole domain. The approach is based on a mixed reformulation of the governing equation and utilizes the standard least-squares objective, with deep neural networks as ansatz functions to approximate the conductivity and flux simultaneously. We provide a thorough analysis of the deep neural network approximations of the conductivity for both continuous and empirical losses, including rigorous error estimates that are explicit in terms of the noise level, various penalty parameters and neural network architectural parameters (depth, width and parameter bound). We also provide multiple numerical experiments in two- and multi-dimensions to illustrate distinct features of the approach, e.g., excellent stability with respect to data noise and capability of solving high-dimensional problems.

Marco Nittscher, Michael Lameter, Riccardo Barbano et al.

The deep image prior (DIP) is a well-established unsupervised deep learning method for image reconstruction; yet it is far from being flawless. The DIP overfits to noise if not early stopped, or optimized via a regularized objective. We build on the regularized fine-tuning of a pretrained DIP, by adopting a novel strategy that restricts the learning to the adaptation of singular values. The proposed SVD-DIP uses ad hoc convolutional layers whose pretrained parameters are decomposed via the singular value decomposition. Optimizing the DIP then solely consists in the fine-tuning of the singular values, while keeping the left and right singular vectors fixed. We thoroughly validate the proposed method on real-measured $μ$CT data of a lotus root as well as two medical datasets (LoDoPaB and Mayo). We report significantly improved stability of the DIP optimization, by overcoming the overfitting to noise.

Bangti Jin, Zehui Zhou, Jun Zou

Invertible neural networks (INNs) represent an important class of deep neural network architectures that have been widely used in several applications. The universal approximation properties of INNs have also been established recently. However, the approximation rate of INNs is largely missing. In this work, we provide an analysis of the capacity of a class of coupling-based INNs to approximate bi-Lipschitz continuous mappings on a compact domain, and the result shows that it can well approximate both forward and inverse maps simultaneously. Furthermore, we develop an approach for approximating bi-Lipschitz maps on infinite-dimensional spaces that simultaneously approximate the forward and inverse maps, by combining model reduction with principal component analysis and INNs for approximating the reduced map, and we analyze the overall approximation error of the approach. Preliminary numerical results show the feasibility of the approach for approximating the solution operator for parameterized second-order elliptic problems.

Riccardo Barbano, Alexander Denker, Hyungjin Chung et al.

Denoising diffusion models have emerged as the go-to generative framework for solving inverse problems in imaging. A critical concern regarding these models is their performance on out-of-distribution tasks, which remains an under-explored challenge. Using a diffusion model on an out-of-distribution dataset, realistic reconstructions can be generated, but with hallucinating image features that are uniquely present in the training dataset. To address this discrepancy during train-test time and improve reconstruction accuracy, we introduce a novel sampling framework called Steerable Conditional Diffusion. Specifically, this framework adapts the diffusion model, concurrently with image reconstruction, based solely on the information provided by the available measurement. Utilising our proposed method, we achieve substantial enhancements in out-of-distribution performance across diverse imaging modalities, advancing the robust deployment of denoising diffusion models in real-world applications.

Tianhao Hu, Bangti Jin, Zhi Zhou

The inverse problem of recovering point sources represents an important class of applied inverse problems. However, there is still a lack of neural network-based methods for point source identification, mainly due to the inherent solution singularity. In this work, we develop a novel algorithm to identify point sources, utilizing a neural network combined with a singularity enrichment technique. We employ the fundamental solution and neural networks to represent the singular and regular parts, respectively, and then minimize an empirical loss involving the intensities and locations of the unknown point sources, as well as the parameters of the neural network. Moreover, by combining the conditional stability argument of the inverse problem with the generalization error of the empirical loss, we conduct a rigorous error analysis of the algorithm. We demonstrate the effectiveness of the method with several challenging experiments.

Bangti Jin, Zehui Zhou

Stochastic variance reduced gradient (SVRG) is an accelerated version of stochastic gradient descent based on variance reduction, and is promising for solving large-scale inverse problems. In this work, we analyze SVRG and a regularized version that incorporates a priori knowledge of the problem, for solving linear inverse problems in Hilbert spaces. We prove that, with suitable constant step size schedules and regularity conditions, the regularized SVRG can achieve optimal convergence rates in terms of the noise level without any early stopping rules, provided that the truncation level is chosen suitably, and standard SVRG is also optimal for problems with nonsmooth solutions under a priori stopping rules. The analysis is based on an explicit error recursion and suitable a priori estimates on the inner loop updates with respect to the anchor point. Numerical experiments are provided to complement the theoretical analysis.

Derick Nganyu Tanyu, Jianfeng Ning, Andreas Hauptmann et al.

Electrical Impedance Tomography (EIT) is a powerful imaging technique with diverse applications, e.g., medical diagnosis, industrial monitoring, and environmental studies. The EIT inverse problem is about inferring the internal conductivity distribution of an object from measurements taken on its boundary. It is severely ill-posed, necessitating advanced computational methods for accurate image reconstructions. Recent years have witnessed significant progress, driven by innovations in analytic-based approaches and deep learning. This review explores techniques for solving the EIT inverse problem, focusing on the interplay between contemporary deep learning-based strategies and classical analytic-based methods. Four state-of-the-art deep learning algorithms are rigorously examined, harnessing the representational capabilities of deep neural networks to reconstruct intricate conductivity distributions. In parallel, two analytic-based methods, rooted in mathematical formulations and regularisation techniques, are dissected for their strengths and limitations. These methodologies are evaluated through various numerical experiments, encompassing diverse scenarios that reflect real-world complexities. A suite of performance metrics is employed to assess the efficacy of these methods. These metrics collectively provide a nuanced understanding of the methods' ability to capture essential features and delineate complex conductivity patterns. One novel feature of the study is the incorporation of variable conductivity scenarios, introducing a level of heterogeneity that mimics textured inclusions. This departure from uniform conductivity assumptions mimics realistic scenarios where tissues or materials exhibit spatially varying electrical properties. Exploring how each method responds to such variable conductivity scenarios opens avenues for understanding their robustness and adaptability.

Javier Antorán, Riccardo Barbano, Johannes Leuschner et al.

Existing deep-learning based tomographic image reconstruction methods do not provide accurate estimates of reconstruction uncertainty, hindering their real-world deployment. This paper develops a method, termed as the linearised deep image prior (DIP), to estimate the uncertainty associated with reconstructions produced by the DIP with total variation regularisation (TV). Specifically, we endow the DIP with conjugate Gaussian-linear model type error-bars computed from a local linearisation of the neural network around its optimised parameters. To preserve conjugacy, we approximate the TV regulariser with a Gaussian surrogate. This approach provides pixel-wise uncertainty estimates and a marginal likelihood objective for hyperparameter optimisation. We demonstrate the method on synthetic data and real-measured high-resolution 2D $μ$CT data, and show that it provides superior calibration of uncertainty estimates relative to previous probabilistic formulations of the DIP. Our code is available at https://github.com/educating-dip/bayes_dip.

Bangti Jin, Jing Li, Yifeng Xu et al.

In this work, we develop an adaptive nonconforming finite element algorithm for the numerical approximation of phase-field parameterized topology optimization governed by the Stokes system. We employ the conforming linear finite element space to approximate the phase field, and the nonconforming linear finite elements (Crouzeix-Raviart elements) and piecewise constants to approximate the velocity field and the pressure field, respectively. We establish the convergence of the adaptive method, i.e., the sequence of minimizers contains a subsequence that converges to a solution of the first-order optimality system, and the associated subsequence of discrete pressure fields also converges. The analysis relies crucially on a new discrete compactness result of nonconforming linear finite elements over a sequence of adaptively generated meshes. We present numerical results for several examples to illustrate the performance of the algorithm, including a comparison with the uniform refinement strategy.

Bangti Jin, Longjun Wu

Physics informed neural networks (PINNs) represent a very popular class of neural solvers for partial differential equations. In practice, one often employs stochastic gradient descent type algorithms to train the neural network. Therefore, the convergence guarantee of stochastic gradient descent is of fundamental importance. In this work, we establish the linear convergence of stochastic gradient descent / flow in training over-parameterized two layer PINNs for a general class of activation functions in the sense of high probability. These results extend the existing result [18] in which gradient descent was analyzed. The challenge of the analysis lies in handling the dynamic randomness introduced by stochastic optimization methods. The key of the analysis lies in ensuring the positive definiteness of suitable Gram matrices during the training. The analysis sheds insight into the dynamics of the optimization process, and provides guarantees on the neural networks trained by stochastic algorithms.

Ruchi Guo, Jiahua Jiang, Bangti Jin et al.

Fluorescence Molecular Tomography (FMT) is a widely used non-invasive optical imaging technology in biomedical research. It usually faces significant accuracy challenges in depth reconstruction, and conventional iterative methods struggle with poor $z$-resolution even with advanced regularization. Supervised learning approaches can improve recovery accuracy but rely on large, high-quality paired training dataset that is often impractical to acquire in practice. This naturally raises the question of how learning-based approaches can be effectively combined with iterative schemes to yield more accurate and stable algorithms. In this work, we present a novel warm-basis iterative projection method (WB-IPM) and establish its theoretical underpinnings. The method is able to achieve significantly more accurate reconstructions than the learning-based and iterative-based methods. In addition, it allows a weaker loss function depending solely on the directional component of the difference between ground truth and neural network output, thereby substantially reducing the training effort. These features are justified by our error analysis as well as simulated and real-data experiments.

Tianhao Hu, Bangti Jin, Fengru Wang

In this work, we present a novel iterative deep Ritz method (IDRM) for solving a general class of elliptic problems. It is inspired by the iterative procedure for minimizing the loss during the training of the neural network, but at each step encodes the geometry of the underlying function space and incorporates a convex penalty to enhance the performance of the algorithm. The algorithm is applicable to elliptic problems involving a monotone operator (not necessarily of variational form) and does not impose any stringent regularity assumption on the solution. It improves several existing neural PDE solvers, e.g., physics informed neural network and deep Ritz method, in terms of the accuracy for the concerned class of elliptic problems. Further, we establish a convergence rate for the method using tools from geometry of Banach spaces and theory of monotone operators, and also analyze the learning error. To illustrate the effectiveness of the method, we present several challenging examples, including a comparative study with existing techniques.

Riccardo Barbano, Johannes Leuschner, Maximilian Schmidt et al.

Deep image prior (DIP) was recently introduced as an effective unsupervised approach for image restoration tasks. DIP represents the image to be recovered as the output of a deep convolutional neural network, and learns the network's parameters such that the output matches the corrupted observation. Despite its impressive reconstructive properties, the approach is slow when compared to supervisedly learned, or traditional reconstruction techniques. To address the computational challenge, we bestow DIP with a two-stage learning paradigm: (i) perform a supervised pretraining of the network on a simulated dataset; (ii) fine-tune the network's parameters to adapt to the target reconstruction task. We provide a thorough empirical analysis to shed insights into the impacts of pretraining in the context of image reconstruction. We showcase that pretraining considerably speeds up and stabilizes the subsequent reconstruction task from real-measured 2D and 3D micro computed tomography data of biological specimens. The code and additional experimental materials are available at https://educateddip.github.io/docs.educated_deep_image_prior/.

Chen Zhang, Riccardo Barbano, Bangti Jin

Learned image reconstruction techniques using deep neural networks have recently gained popularity, and have delivered promising empirical results. However, most approaches focus on one single recovery for each observation, and thus neglect the uncertainty information. In this work, we develop a novel computational framework that approximates the posterior distribution of the unknown image at each query observation. The proposed framework is very flexible: It handles implicit noise models and priors, it incorporates the data formation process (i.e., the forward operator), and the learned reconstructive properties are transferable between different datasets. Once the network is trained using the conditional variational autoencoder loss, it provides a computationally efficient sampler for the approximate posterior distribution via feed-forward propagation, and the summarizing statistics of the generated samples are used for both point-estimation and uncertainty quantification. We illustrate the proposed framework with extensive numerical experiments on positron emission tomography (with both moderate and low count levels) showing that the framework generates high-quality samples when compared with state-of-the-art methods.

Riccardo Barbano, Zeljko Kereta, Andreas Hauptmann et al.

Deep learning-based image reconstruction approaches have demonstrated impressive empirical performance in many imaging modalities. These approaches usually require a large amount of high-quality paired training data, which is often not available in medical imaging. To circumvent this issue we develop a novel unsupervised knowledge-transfer paradigm for learned reconstruction within a Bayesian framework. The proposed approach learns a reconstruction network in two phases. The first phase trains a reconstruction network with a set of ordered pairs comprising of ground truth images of ellipses and the corresponding simulated measurement data. The second phase fine-tunes the pretrained network to more realistic measurement data without supervision. By construction, the framework is capable of delivering predictive uncertainty information over the reconstructed image. We present extensive experimental results on low-dose and sparse-view computed tomography showing that the approach is competitive with several state-of-the-art supervised and unsupervised reconstruction techniques. Moreover, for test data distributed differently from the training data, the proposed framework can significantly improve reconstruction quality not only visually, but also quantitatively in terms of PSNR and SSIM, when compared with learned methods trained on the synthetic dataset only.

Riccardo Barbano, Željko Kereta, Chen Zhang et al.

Image reconstruction methods based on deep neural networks have shown outstanding performance, equalling or exceeding the state-of-the-art results of conventional approaches, but often do not provide uncertainty information about the reconstruction. In this work we propose a scalable and efficient framework to simultaneously quantify aleatoric and epistemic uncertainties in learned iterative image reconstruction. We build on a Bayesian deep gradient descent method for quantifying epistemic uncertainty, and incorporate the heteroscedastic variance of the noise to account for the aleatoric uncertainty. We show that our method exhibits competitive performance against conventional benchmarks for computed tomography with both sparse view and limited angle data. The estimated uncertainty captures the variability in the reconstructions, caused by the restricted measurement model, and by missing information, due to the limited angle geometry.

Riccardo Barbano, Chen Zhang, Simon Arridge et al.

Recent advances in reconstruction methods for inverse problems leverage powerful data-driven models, e.g., deep neural networks. These techniques have demonstrated state-of-the-art performances for several imaging tasks, but they often do not provide uncertainty on the obtained reconstruction. In this work, we develop a scalable, data-driven, knowledge-aided computational framework to quantify the model uncertainty via Bayesian neural networks. The approach builds on, and extends deep gradient descent, a recently developed greedy iterative training scheme, and recasts it within a probabilistic framework. Scalability is achieved by being hybrid in the architecture: only the last layer of each block is Bayesian, while the others remain deterministic, and by being greedy in training. The framework is showcased on one representative medical imaging modality, viz. computed tomography with either sparse view or limited view data, and exhibits competitive performance with respect to state-of-the-art benchmarks, e.g., total variation, deep gradient descent and learned primal-dual.

Chen Zhang, Bangti Jin

Aleatoric uncertainty is an intrinsic property of ill-posed inverse and imaging problems. Its quantification is vital for assessing the reliability of relevant point estimates. In this paper, we propose an efficient framework for quantifying aleatoric uncertainty for deep residual learning and showcase its significant potential on image restoration. In the framework, we divide the conditional probability modeling for the residual variable into a deterministic homo-dimensional level, a stochastic low-dimensional level and a merging level. The low-dimensionality is especially suitable for sparse correlation between image pixels, enables efficient sampling for high dimensional problems and acts as a regularizer for the distribution. Preliminary numerical experiments show that the proposed method can give not only state-of-the-art point estimates of image restoration but also useful associated uncertainty information.

Bangti Jin, Yubin Yan, Zhi Zhou

We develop and analyze a numerical method for stochastic time-fractional diffusion driven by additive fractionally integrated Gaussian noise. The model involves two nonlocal terms in time, i.e., a Caputo fractional derivative of order $α\in(0,1)$, and fractionally integrated Gaussian noise (with a Riemann-Liouville fractional integral of order $γ\in[0,1]$ in the front). The numerical scheme approximates the model in space by the Galerkin method with continuous piecewise linear finite elements and in time by the classical Grünwald-Letnikov method, and the noise by the $L^2$-projection. Sharp strong and weak convergence rates are established, using suitable nonsmooth data error estimates for the deterministic counterpart. Numerical results are presented to support the theoretical findings.

Bangti Jin, Buyang Li, Zhi Zhou

In this work, a complete error analysis is presented for fully discrete solutions of the subdiffusion equation with a time-dependent diffusion coefficient, obtained by the Galerkin finite element method with conforming piecewise linear finite elements in space and backward Euler convolution quadrature in time. The regularity of the solutions of the subdiffusion model is proved for both nonsmooth initial data and incompatible source term. Optimal-order convergence of the numerical solutions is established using the proven solution regularity and a novel perturbation argument via freezing the diffusion coefficient at a fixed time. The analysis is supported by numerical experiments.

Beiping Duan, Bangti Jin, Raytcho Lazarov et al.

We present and analyze a space-time Petrov-Galerkin finite element method for a time-fractional diffusion equation involving a Riemann-Liouville fractional derivative of order $α\in(0,1)$ in time and zero initial data. We derive a proper weak formulation involving different solution and test spaces and show the inf-sup condition for the bilinear form and thus its well-posedness. Further, we develop a novel finite element formulation, show the well-posedness of the discrete problem, and establish error bounds in both energy and $L^2$ norms for the finite element solution. In the proof of the discrete inf-sup condition, a certain nonstandard $L^2$ stability property of the $L^2$ projection operator plays a key role. We provide extensive numerical examples to verify the convergence of the method.

Eric T. Chung, Yalchin Efendiev, Bangti Jin et al.

In this work, we propose a generalized multiscale inversion algorithm for heterogeneous problems that aims at solving an inverse problem on a computational coarse grid. Previous inversion techniques for multiscale problems seek a coarse-grid media properties, e.g., permeability and conductivity, and by doing so, they assume that there exists a homogenized representation of the underlying fine-scale permeability field on a coarse grid. Generally such assumptions do not hold for highly heterogeneous fields, e.g., fracture media or channelized fields, where the width of channels are very small compared to the coarse-grid sizes. In these cases, grid refinement can lead to many degrees of freedom, and thus unattractive to apply. The proposed algorithm is based on the Generalized Multiscale Finite Element Method (GMsFEM), which uses local spectral problems to identify non-localized features, i.e., channels (high-conductivity inclusions that connect the boundaries of the coarse-grid block). The inclusion of these features in the coarse space enables one to achieve a good accuracy. The approach is valid under the assumption that the solution can be well represented in a reduced-dimensional space by multiscale basis functions. In practice, these basis functions are non-observable as we do not identify the fine-scale features of the permeability field. Our inversion algorithm finds the discretization parameters of the resulting system. By doing so, we identify the appropriate coarse-grid parameters representing the permeability field instead of fine-grid permeability field. We illustrate the approach by numerical results for fractured media.

Bangti Jin, Raytcho Lazarov, Zhi Zhou

We consider the initial/boundary value problem for the fractional diffusion and diffusion-wave equations involving a Caputo fractional derivative in time. We develop two "simple" fully discrete schemes based on the Galerkin finite element method in space and convolution quadrature in time with the generating function given by the implicit backward Euler method/second-order backward difference method, and establish error estimates optimal with respect to the regularity of the initial data. These two schemes are first and second-order accurate in time for nonsmooth initial data. Extensive numerical experiments for one and two-dimensional problems confirm the convergence analysis. A detailed comparison with several popular time stepping schemes is also performed. The numerical results indicate that the proposed fully discrete schemes are accurate and robust for nonsmooth data, and competitive with existing schemes.

Bangti Jin, Raytcho Lazarov, Vidar Thomée et al.

We consider three types of subdiffusion models, namely single-term, multi-term and distributed order fractional diffusion equations, for which the maximum-principle holds and which, in particular, preserve nonnegativity. Hence the solution is nonnegative for nonnegative initial data. Following earlier work on the heat equation, our purpose is to study whether this property is inherited by certain spatially semidiscrete and fully discrete piecewise linear finite element methods, including the standard Galerkin method, the lumped mass method and the finite volume element method. It is shown that, as for the heat equation, when the mass matrix is nondiagonal, nonnegativity is not preserved for small time or time-step, but may reappear after a positivity threshold. For the lumped mass method nonnegativity is preserved if and only if the triangulation in the finite element space is of Delaunay type. Numerical experiments illustrate and complement the theoretical results.

Bangti Jin, Raytcho Lazarov, Dongwoo Sheen et al.

In this work, we consider the numerical solution of an initial boundary value problem for the distributed order time fractional diffusion equation. The model arises in the mathematical modeling of ultra-slow diffusion processes observed in some physical problems, whose solution decays only logarithmically as the time $t$ tends to infinity. We develop a space semidiscrete scheme based on the standard Galerkin finite element method, and establish error estimates optimal with respect to data regularity in $L^2(D)$ and $H^1(D)$ norms for both smooth and nonsmooth initial data. Further, we propose two fully discrete schemes, based on the Laplace transform and convolution quadrature generated by the backward Euler method, respectively, and provide optimal convergence rates in the $L^2(D)$ norm, which exhibits exponential convergence and first-order convergence in time, respectively. Extensive numerical experiments are provided to verify the error estimates for both smooth and nonsmooth initial data, and to examine the asymptotic behavior of the solution.

Yuling Jiao, Bangti Jin, Xiliang Lu

We develop a primal dual active set with continuation algorithm for solving the \ell^0-regularized least-squares problem that frequently arises in compressed sensing. The algorithm couples the the primal dual active set method with a continuation strategy on the regularization parameter. At each inner iteration, it first identifies the active set from both primal and dual variables, and then updates the primal variable by solving a (typically small) least-squares problem defined on the active set, from which the dual variable can be updated explicitly. Under certain conditions on the sensing matrix, i.e., mutual incoherence property or restricted isometry property, and the noise level, the finite step global convergence of the algorithm is established. Extensive numerical examples are presented to illustrate the efficiency and accuracy of the algorithm and the convergence analysis.

Jian Huang, Yuling Jiao, Bangti Jin et al.

In this paper, we consider the problem of recovering a sparse signal based on penalized least squares formulations. We develop a novel algorithm of primal-dual active set type for a class of nonconvex sparsity-promoting penalties, including $\ell^0$, bridge, smoothly clipped absolute deviation, capped $\ell^1$ and minimax concavity penalty. First we establish the existence of a global minimizer for the related optimization problems. Then we derive a novel necessary optimality condition for the global minimizer using the associated thresholding operator. The solutions to the optimality system are coordinate-wise minimizers, and under minor conditions, they are also local minimizers. Upon introducing the dual variable, the active set can be determined using the primal and dual variables together. Further, this relation lends itself to an iterative algorithm of active set type which at each step involves first updating the primal variable only on the active set and then updating the dual variable explicitly. When combined with a continuation strategy on the regularization parameter, the primal dual active set method is shown to converge globally to the underlying regression target under certain regularity conditions. Extensive numerical experiments with both simulated and real data demonstrate its superior performance in efficiency and accuracy compared with the existing sparse recovery methods.