If I recall correctly, Matlab's GMM regularization works by adding a multiple of the identity matrix to the covariance matrices.
Using regularization can be viewed as placing a prior on the covariance matrices. As the strength of regularization increases, the model becomes less flexible. As a consequence, the training set log likelihood may be lower for more strongly regularized models. Because of this, training set log likelihood isn't a great way to compare models with different levels of regularization.
AIC and BIC penalize more flexible models, so this is a step in the right direction. They do this by measuring the number of free parameters, and are straightforward to use for comparing models with different number of mixture components and different types of covariance matrices (tied, diagonal, full, etc.). In these cases, the number of free parameters can simply be counted. The situation is more complicated for regularized covariance matrices, because model flexibility is governed by a continuous regularization parameter, rather than a change in the number of parameters. It may be possible to find an expression for the 'effective degrees of freedom' as a function of the regularization parameter. For example, this can be done in the case of ridge regression.
Two alternatives come to mind (I'm assuming you're holding the number of components, type of covariance matrix, etc. fixed, and only changing the regularization parameter). 1) Train each model on a subset of your data, then compute the log likelihood on the remaining test set (better yet, use cross validation). Choose the regularization parameter that maximizes test set likelihood. 2) Use Bayesian model comparison. Choose the regularization parameter that maximizes the evidence. This is a more complicated approach, because you'd have to define priors on all parameters, derive an expression for the evidence, and find a way to compute or approximate it.
One reason to prefer cross validation/Bayesian model comparison over AIC/BIC is that AIC/BIC rely on asymptotic approximations, which may or may not be appropriate for your data. Note that cross validation is not a valid strategy for choosing the number of mixture components, because the test set likelihood may increase monotonically with the number of components.