A scalable algorithm for MAP estimators in Bayesian inverse problems with Besov priors

We present a scalable solver for approximating the maximum a posteriori (MAP) point of Bayesian inverse problems with Besov priors based on wavelet expansions with random coefficients. It is a subspace trust region interior reflective Newton conjugate gradient method for bound constrained optimization problems. The method combines the rapid locally-quadratic convergence rate properties of Newton's method, the effectiveness of trust region globalization for treating ill-conditioned problems, and the Eisenstat--Walker idea of preventing oversolving. We demonstrate the scalability of the proposed method on two inverse problems: a deconvolution problem and a coefficient inverse problem governed by elliptic partial differential equations. The numerical results show that the number of Newton iterations is independent of the number of wavelet coefficients $n$ and the computation time scales linearly in $n$. It will be numerically shown, under our implementations, that the proposed solver is two times faster than the split Bregman approach, and it is an order of magnitude less expensive than the interior path following primal-dual method. Our results also confirm the fact that the Besov $\mathbb{B}_{11}^1$ prior is sparsity promoting, discretization-invariant, and edge-preserving for both imaging and inverse problems governed by partial differential equations.