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 ofBayesian inverse problems with Besov priors based on waveletexpansions with random coefficients. It is a subspace trust regioninterior reflective Newton conjugate gradient method for boundconstrained optimization problems. The method combines the rapidlocally-quadratic convergence rate properties of Newton's method, theeffectiveness of trust region globalization for treatingill-conditioned problems, and the Eisenstat--Walker idea of preventingoversolving. We demonstrate the scalability of the proposed method ontwo inverse problems: a deconvolution problem and a coefficientinverse problem governed by elliptic partial differentialequations. The numerical results show that the number of Newtoniterations is independent of the number of wavelet coefficients nand 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 splitBregman approach, and it is an order of magnitude less expensive thanthe interior path following primal-dual method. Our results also confirm the fact that the Besov B₁₁¹ prior is sparsity promoting, discretization-invariant, and edge-preserving for both imaging and inverse problems governed by partial differential equations.