Numerical Finite-element method (FEM) based algorithms have been widely applied for the reconstruction of the photoacoustic image. As compared with the traditional analytic methods, the FEM based methods can be easily used to deal with problems with irregularly shaped imaging domain. However, the FEM based algorithms are usually computationally intensive because repeated manipulations of matrices with larger size are needed during the reconstruction process. To tackle such a problem, a novel method is proposed for reducing the size of the matrix to be inversed during the reconstruction process and hence speed up the inverse reconstruction without any sacrifice of the reconstruction accuracy.