Computing the regularized solution of Bayesian linear inverse problems as well as the corresponding regularization parameter is highly desirable in many applications. This paper proposes a novel iterative method, termed the Projected Newton method (PNT), that can simultaneously update the regularization parameter and solution step by step without requiring any high-cost matrix inversions or decompositions. By reformulating the Tikhonov regularization as a constrained minimization problem and writing its Lagrangian function, a Newton-type method coupled with a Krylov subspace method, called the generalized Golub-Kahan bidiagonalization, is employed for the unconstrained Lagrangian function. The resulting PNT algorithm only needs solving a small-scale linear system to get a descent direction of a merit function at each iteration, thus significantly reducing computational overhead. Rigorous convergence results are proved, showing that PNT always converges to the unique regularized solution and the corresponding Lagrangian multiplier. Experimental results on both small and large-scale Bayesian inverse problems demonstrate its excellent convergence property, robustness and efficiency. Given that the most demanding computational tasks in PNT are primarily matrix-vector products, it is particularly well-suited for large-scale problems.
One of the key elements of probabilistic seismic risk assessment studies is the fragility curve, which represents the conditional probability of failure of a mechanical structure for a given scalar measure derived from seismic ground motion. For many structures of interest, estimating these curves is a daunting task because of the limited amount of data available; data which is only binary in our framework, i.e., only describing the structure as being in a failure or non-failure state. A large number of methods described in the literature tackle this challenging framework through parametric log-normal models. Bayesian approaches, on the other hand, allow model parameters to be learned more efficiently. However, the impact of the choice of the prior distribution on the posterior distribution cannot be readily neglected and, consequently, neither can its impact on any resulting estimation. This paper proposes a comprehensive study of this parametric Bayesian estimation problem for limited and binary data. Using the reference prior theory as a cornerstone, this study develops an objective approach to choosing the prior. This approach leads to the Jeffreys prior, which is derived for this problem for the first time. The posterior distribution is proven to be proper with the Jeffreys prior but improper with some traditional priors found in the literature. With the Jeffreys prior, the posterior distribution is also shown to vanish at the boundaries of the parameters' domain, which means that sampling the posterior distribution of the parameters does not result in anomalously small or large values. Therefore, the use of the Jeffreys prior does not result in degenerate fragility curves such as unit-step functions, and leads to more robust credibility intervals. The numerical results obtained from different case studies-including an industrial example-illustrate the theoretical predictions.
A bottleneck of sufficient dimension reduction (SDR) in the modern era is that, among numerous methods, only the sliced inverse regression (SIR) is generally applicable under the high-dimensional settings. The higher-order inverse regression methods, which form a major family of SDR methods that are superior to SIR in the population level, suffer from the dimensionality of their intermediate matrix-valued parameters that have an excessive number of columns. In this paper, we propose the generic idea of using a small subset of columns of the matrix-valued parameter for SDR estimation, which breaks the convention of using the ambient matrix for the higher-order inverse regression methods. With the aid of a quick column selection procedure, we then generalize these methods as well as their ensembles towards sparsity under the ultrahigh-dimensional settings, in a uniform manner that resembles sparse SIR and without additional assumptions. This is the first promising attempt in the literature to free the higher-order inverse regression methods from their dimensionality, which facilitates the applicability of SDR. The gain of column selection with respect to SDR estimation efficiency is also studied under the fixed-dimensional settings. Simulation studies and a real data example are provided at the end.
The method of fundamental solutions (MFS), also known as the method of auxiliary sources (MAS), is a well-known computational method for the solution of boundary-value problems. The final solution ("MAS solution") is obtained once we have found the amplitudes of $N$ auxiliary "MAS sources." Past studies have demonstrated that it is possible for the MAS solution to converge to the true solution even when the $N$ auxiliary sources diverge and oscillate. The present paper extends the past studies by demonstrating this possibility within the context of Laplace's equation with Neumann boundary conditions. One can thus obtain the correct solution from sources that, when $N$ is large, must be considered unphysical. We carefully explain the underlying reasons for the unphysical results, distinguish from other difficulties that might concurrently arise, and point to significant differences with time-dependent problems that were studied in the past.
When an inverse problem is solved by a gradient-based optimization algorithm, the corresponding forward and adjoint problems, which are introduced to compute the gradient, can be also solved iteratively. The idea of iterating at the same time on the inverse problem unknown and on the forward and adjoint problem solutions yields the concept of one-shot inversion methods. We are especially interested in the case where the inner iterations for the direct and adjoint problems are incomplete, that is, stopped before achieving a high accuracy on their solutions. Here, we focus on general linear inverse problems and generic fixed-point iterations for the associated forward problem. We analyze variants of the so-called multi-step one-shot methods, in particular semi-implicit schemes with a regularization parameter. We establish sufficient conditions on the descent step for convergence, by studying the eigenvalues of the block matrix of the coupled iterations. Several numerical experiments are provided to illustrate the convergence of these methods in comparison with the classical gradient descent, where the forward and adjoint problems are solved exactly by a direct solver instead. We observe that very few inner iterations are enough to guarantee good convergence of the inversion algorithm, even in the presence of noisy data.
Sparse recovery principles play an important role in solving many nonlinear ill-posed inverse problems. We investigate a variational framework with support Oracle for compressed sensing sparse reconstructions, where the available measurements are nonlinear and possibly corrupted by noise. A graph neural network, named Oracle-Net, is proposed to predict the support from the nonlinear measurements and is integrated into a regularized recovery model to enforce sparsity. The derived nonsmooth optimization problem is then efficiently solved through a constrained proximal gradient method. Error bounds on the approximate solution of the proposed Oracle-based optimization are provided in the context of the ill-posed Electrical Impedance Tomography problem. Numerical solutions of the EIT nonlinear inverse reconstruction problem confirm the potential of the proposed method which improves the reconstruction quality from undersampled measurements, under sparsity assumptions.
A moving mesh finite element method is studied for the numerical solution of Bernoulli free boundary problems. The method is based on the pseudo-transient continuation with which a moving boundary problem is constructed and its steady-state solution is taken as the solution of the underlying Bernoulli free boundary problem. The moving boundary problem is solved in a split manner at each time step: the moving boundary is updated with the Euler scheme, the interior mesh points are moved using a moving mesh method, and the corresponding initial-boundary value problem is solved using the linear finite element method. The method can take full advantages of both the pseudo-transient continuation and the moving mesh method. Particularly, it is able to move the mesh, free of tangling, to fit the varying domain for a variety of geometries no matter if they are convex or concave. Moreover, it is convergent towards steady state for a broad class of free boundary problems and initial guesses of the free boundary. Numerical examples for Bernoulli free boundary problems with constant and non-constant Bernoulli conditions and for nonlinear free boundary problems are presented to demonstrate the accuracy and robustness of the method and its ability to deal with various geometries and nonlinearities.
We present and analyze an a posteriori error estimator for a space-time hybridizable discontinuous Galerkin discretization of the time-dependent advection-diffusion problem. The residual-based error estimator is proven to be reliable and locally efficient. In the reliability analysis we combine a Peclet-robust coercivity type result and a saturation assumption, while local efficiency analysis is based on using bubble functions. The analysis considers both local space and time adaptivity and is verified by numerical simulations on problems which include boundary and interior layers.
We study a truncated two-dimensional moment problem in terms of the Stieltjes transform. The set of the solutions is described by the Schur step-by-step algorithm, which is based on the continued fraction expansion of the solution. In particular, the obtained results are applicable to the two-dimensional moment problem for atomic measures.
In this paper we consider ill-posed inverse problems, both linear and nonlinear, by a heavy ball method in which a strongly convex regularization function is incorporated to detect the feature of the sought solution. We develop ideas on how to adaptively choose the step-sizes and the momentum coefficients to achieve acceleration over the Landweber-type method. We then analyze the method and establish its regularization property when it is terminated by the discrepancy principle. Various numerical results are reported which demonstrate the superior performance of our method over the Landweber-type method by reducing substantially the required number of iterations and the computational time.
An interesting case of the well-known Dataset Shift Problem is the classification of Electroencephalogram (EEG) signals in the context of Brain-Computer Interface (BCI). The non-stationarity of EEG signals can lead to poor generalisation performance in BCI classification systems used in different sessions, also from the same subject. In this paper, we start from the hypothesis that the Dataset Shift problem can be alleviated by exploiting suitable eXplainable Artificial Intelligence (XAI) methods to locate and transform the relevant characteristics of the input for the goal of classification. In particular, we focus on an experimental analysis of explanations produced by several XAI methods on an ML system trained on a typical EEG dataset for emotion recognition. Results show that many relevant components found by XAI methods are shared across the sessions and can be used to build a system able to generalise better. However, relevant components of the input signal also appear to be highly dependent on the input itself.