Due to multi-exponential decay-time properties of the subsurface volume units or layers, magnetic resonance sounding (MRS) relaxation data exhibit a multi-exponential behavior. MRS inverse problem in a multiexponential modeling framework brings about a very large size of the parameter space which is computationally costly. In this paper, a fast and memory efficient inversion algorithm to retrieve the aquifer properties in terms of water content and relaxation time is presented. The original nonsymmetric linearized forward matrix is projected onto a Krylov subspace with smaller dimension using an iterative Golub-Kahan-Lanczos bidiagonalization (GKL) method. Because of ill-conditioning of the projected linearized forward matrix a regularized damped least squares equation is applied at each step of the GKL factorization method to extract the best possible approximation of the partial water content. Numerical experiments based on synthetic and field data demonstrate that the proposed inversion method provides a good estimation of the water content and relaxation time compared to the standard algorithm with computationally more efficient functionality. (C) 2020 Elsevier B.V. All rights reserved.