Functional near infra-red spectroscopy (fNIRs) with near infra-red light have been active research areas for both clinical and pre-clinical applications for more than three decades. The development of more advanced image reconstruction methods is required to improve the accuracy fNIRs of complex tissue structures. In this paper, the least square QR decomposition (LSQR) method for solving the inverse problem has been implemented for real fNIRs data based on working memory (WM). The sensitivity matrix is being generated using the Monte Carlo (MC) simulation. For image reconstruction, the numerical algorithm for the LSQR method is created and implemented in MATLAB. Lastly, the variation of oxy and deoxy haemoglobin levels is monitored based on absorption changes, and the findings obtained using the LSQR regularization method are in good agreement with the real fNIRs WM data. © 2021 IEEE.