  We introduce a distributionally robust maximum likelihood estimation model
with a Wasserstein ambiguity set to infer the inverse covariance matrix of a
$p$-dimensional Gaussian random vector from $n$ independent samples. The
proposed model minimizes the worst case (maximum) of Stein's loss across all
normal reference distributions within a prescribed Wasserstein distance from
the normal distribution characterized by the sample mean and the sample
covariance matrix. We prove that this estimation problem is equivalent to a
semidefinite program that is tractable in theory but beyond the reach of
general purpose solvers for practically relevant problem dimensions $p$. In the
absence of any prior structural information, the estimation problem has an
analytical solution that is naturally interpreted as a nonlinear shrinkage
estimator. Besides being invertible and well-conditioned even for $p>n$, the
new shrinkage estimator is rotation-equivariant and preserves the order of the
eigenvalues of the sample covariance matrix. These desirable properties are not
imposed ad hoc but emerge naturally from the underlying distributionally robust
optimization model. Finally, we develop a sequential quadratic approximation
algorithm for efficiently solving the general estimation problem subject to
conditional independence constraints typically encountered in Gaussian
graphical models.
