We have developed a 2D anisotropic magnetotelluric (MT) inversion algorithm that uses a limited-memory quasi-Newton (QN) method for bounds-constrained optimization. This algorithm solves the inverse problem, which is nonlinear, by iterative minimization of linearized approximations of the classic Tikhonov regularized objective function. The QN approximation for the Hessian matrix is only implemented for the data-misfit term of the objective function; the part of the Hessian matrix for the regularization is explicitly computed. This adjustment results in a better approximation for the data-misfit term in particular. The inversion algorithm considers arbitrary anisotropy, and it is extended for special cases including azimuthal and vertical anisotropy. The algorithm is shown to be stable and converges rapidly for several simple anisotropic models. These synthetic tests also confirm that the anisotropic inversion produces a correct anomaly with different but equivalent anisotropic parameters. We also consider a complex 2D anisotropic model; the successful results for this model further confirm that the inversion algorithm presented here, which uses the novel modified limited-memory QN approach, is capable of solving the 2D anisotropic MT inverse problem. Finally, we have evaluated a practical application on MT data collected in northern Tibet to demonstrate the effectiveness and stability of our algorithm.