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 the weight matrix are all represented by wavelet expansions, thus yielding a multiresolution 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's) 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 the conjugate gradient descent method. This approach has been applied to continuous wave data calculated based on the diffusion approximation of several two-dimensional (2-D) test media. Compared to a previously reported one grid algorithm, the multigrid method requires substantially shorter computation time under the same reconstruction quality criterion.