Unfolding problems often arise in the context of statistical data analysis. Such problematics occur when the probability distribution of a physical quantity is to be measured, but it is randomized (smeared) by some well understood process, such as a non-ideal detector response or a well described physical phenomenon. In such case it is said that the original probability distribution of interest is folded by a known response function. The reconstruction of the original probability distribution from the measured one is called unfolding. That technically involves evaluation of the non-bounded inverse of an integral operator over the space of L^1 functions, which is known to be an ill-posed problem. For the pertinent regularized operator inversion, we propose a linear iterative formula and provide proof of convergence in a probability theory context. Furthermore, we provide formulae for error estimates at finite iteration stopping order which are of utmost importance in practical applications: the approximation error, the propagated statistical error, and the propagated systematic error can be quantified. The arguments are based on the Riesz-Thorin theorem mapping the original L^1 problem to L^2 space, and subsequent application of ordinary L^2 spectral theory of operators. A library implementation in C of the algorithm along with corresponding error propagation is also provided. A numerical example also illustrates the method in operation.