Aims. In this study, we compare several methods of modeling large-scale systematic differences between catalogs of positions of extragalactic radio sources provided by very long baseline interferometry with an emphasis on mitigating the impact of outliers. Methods. The coordinate difference between catalogs was parameterized by first and second order coefficients of vector spherical harmonics. We solved for these coefficients by using the least-squares method (L2-norm) and, alternatively, by L1-norm minimization. The problem of outliers was addressed either by rejecting them on the basis of their difference to the mean or by using the cell median (CM) method, consisting of reducing the difference field to median values in cells of equal area. The methods were applied to simulated catalogs exhibiting systematics of within 100 microarcseconds - for which we knew the expected results - and to real catalogs. Results. In simulations, the L1 minimization appears practically insensitive to outliers and is within a few microarcseconds of the expected results. Least-squared fitting preceded by L1-norm-based outlier detection performed similarly. The CM method gets close to the true parameter values, within one microarcsecond. When applied to real catalogs, all methods provide close results within a few microarcseconds. Conclusions. The study shows that all tested methods are consistent with each other within a few microarcseconds. Hybrid L2/L1 and iterative L2 methods proved to be very effective in eliminating outliers and showed the best accuracy of the estimated parameters of the mutual orientation of celestial reference frames.