This paper describes an adaptation of the Optimally Localized Averaging (OLA) inversion technique, originally developed for geo- and helioseismological applications, to the interpretation of solar spectroscopic data. It focuses on inverting the thermodynamical properties of the solar atmosphere, assuming that the atmosphere and radiation field are in local thermodynamic equilibrium (LTE). We leave inversions of magnetic field and non-LTE inversions for future work. The advantage with the OLA method is that it computes solutions that are optimally depth resolved with minimal crosstalk error between variables. Additionally, the method allows for direct assessment of the vertical resolution of the inverted solutions. The primary challenges faced when adapting the method to spectroscopic inversions originate with the possible large-amplitude differences between the atmospheric model used to initiate the inversion and the underlying atmosphere it aims to recover, necessitating the development of an iterative scheme. Here, we describe the iterative OLA method we have developed for both single and multivariable inversions and demonstrate its performance on simulated data and synthesized spectra. We note that, when carrying out multivariable inversions, employing response function amplification factors can address the inherent spectral sensitivity bias that makes it hard to invert for less spectrally sensitive variables. The OLA method can, in most cases, reliably invert as well as or better than the frequently employed Stokes Inversion based on Response functions (SIR) scheme, but some difficulties remain. In particular, the method struggles to recover large-scale offsets in the atmospheric stratification. We propose future strategies to improve this aspect.