The model-based image reconstruction techniques for photoacoustic (PA) tomography require an explicit regularization. An error estimate () minimization-based approach was proposed and developed for the determination of a regularization parameter for PA imaging. The regularization was used within Lanczos bidiagonalization framework, which provides the advantage of dimensionality reduction for a large system of equations. It was shown that the proposed method is computationally faster than the state-of-the-art techniques and provides similar performance in terms of quantitative accuracy in reconstructed images. It was also shown that the error estimate () can also be utilized in determining a suitable regularization parameter for other popular techniques such as Tikhonov, exponential, and nonsmooth ( and total variation norm based) regularization methods.