In this paper, we present a wavelet based multigrid approach to solve the perturbation equation encountered in optical tomography. With this scheme, the unknown image, the data, as well as weight matrix are all represented by wavelet expansions, and thus yielding a multi–resoltuion representation of the original perturbation equation in the wavelet domain. This transformed equation is then solved using a multigrid scheme, by which an increasing portion of wavelet coefficients of the unknown image are solved in successive approximations. One can also quickly identify regions of interest (ROI) from a coarse level reconstruction and restrict the reconstruction in the following fine resolutions to those regions. At each resolution level, a regularized least squares solution is obtained using a Conjugate Gradient Descent (CGD) method. Compared to a previously reported one grid algorithm, the multigrid method requires substantially shorter computation time under the same reconstruction quality criterion.