Over the years, seismic amplitude variation with offset has been successfully applied for predicting the elastic properties of the subsurface. Nevertheless, the solution of the amplitude inversion is not stable due to insufficient information in the recorded seismic data. To enhance the stability, various amplitude inversions under a Bayesian framework have been introduced and most of which are based on single-component a priori probability densities, which can be solved deterministically. For multicomponent a priori densities, Monte Carlo sampling methods are often used to obtain the posterior distribution of the inverted parameters. The variational methods, on the other hand, provide the gradients of the parameters with respect to the data. Here, we present an inversion framework that uses the gradients of the variational inference as regularization. Due to the ill-posedness of the amplitude inversion, solving the likelihood of the data leads to solutions that may exist within the low-probability regions of the a priori density. The rule of the regularization (variational inference gradient) is, therefore, to force solutions within the high-probability regions that enhance the stability of the inversion process. The performance of the regularization is tested using synthetic and field data examples to jointly estimate the elastic and Thomsen's anisotropy parameters. The proposed constrained optimization successfully predicts the P-wave velocity, S-wave velocity and P-wave anisotropy and preferably constrains the density and near-vertical anisotropy parameters.