This research is motivated by challenges in addressing optimization models arising in Big Data. Such models are often formulated as large scale stochastic optimization problems. When the probability distribution of the data is unknown, the Sample Average Approximation (SAA) scheme can be employed which results in an Empirical Risk Minimization (EMR) problem. To address this class of problems deterministic solution methods, such as the Broyden, Fletcher, Goldfarb, Shanno (BFGS) method, face high computational cost per iteration and memory requirement issues due to presence of uncertainty and high dimensionality of the solution space. To cope with these challenges, stochastic methods with limited memory variants have been developed recently. However, the solutions generated by such methods might be dense requiring high memory capacity. To generate sparse solutions, in the literature, standard 𝐿1 regularization technique is employed. The regularization term which includes the 𝐿1 regularization parameter and 𝐿1 norm is added to the objective function of the problem. Under this approach, addition of constant 𝐿1 regularization parameter changes the original problem and the solutions obtained by solving the regularized problem are approximate solutions. Moreover, limited information is available in the literature to obtain sparse solutions to the original problem. To address this gap, in this research we develop an iterative 𝐿1 Regularized Limited memory Stochastic BFGS (iRLS-BFGS) method in which the 𝐿1 regularization parameter and the step-size parameter are simultaneously updated at each iteration. Our goal is to find the suitable decay rates for these two sequences in our algorithm. To address this research question, we first implement the iRLS-BFGS algorithm on a Big Data text classification problem and provide a detailed numerical comparison of the performance of the developed algorithm under different choices of the update rules. Our numerical experiments imply that when both the stepsize and the 𝐿1 regularization parameter decay at the rate of the order 1/sqrt(k), the best convergence is achieved. Later, to support our findings, we apply our method to address a large scale image deblurring problem arising in signal processing using the update rule from the previous application. As a result, we obtain much clear deblurred images compared to the classical algorithm’s deblurred output images when both the step-size and the 𝐿1 regularization parameter decay at the rate of the order 1/sqrt(k).