It is a well known property that $E[|X - p|]$ is minimized for $p$ such that $P(X \leq p) = \frac{1}{2}$ (i.e. $p$ is the median value for the random variable $X$).
Now suppose $\mathbf{X} \sim \mathcal{N} (\pmb{\mu}, \Sigma)$. I'm wonder what the expectation $E[|X - p|]$ defined above actually evaluates to. That is, I'm seeking to evaluate $\int_{\infty}^\infty |\mathbf{X} - \mu| f(\mathbf{X}) d\mathbf{x}$. I've played around with the univariate case ($X \sim \mathcal{N} (\mu, \sigma^2)$) a little using some software, and I come up with an interesting closed form expression where the integral comes out to be $\sigma \sqrt{\frac{2}{\pi}}$, but I'm not extremely familiar with Gaussian integrals and I'm not comfortable generalizing this to the multivariate case. Can someone point me in the right direction?
This is a long comment, but it looks like even the $2$-dimensional case will be hard to solve exactly.
Suppose $X$ is $n$-dimensional. Without loss of generality (by translating and rotating $X$), $\mu=0$ and $\Sigma$ is diagonal, giving independent $X_i\sim N(0,\,\sigma_i^2)$. We wish to evaluate$$\int_{\Bbb R^n}\sqrt{\frac{\sum_ix_i^2}{(2\pi)^n\prod_i\sigma_i^2}}\exp\left(-\tfrac12\sum_i\sigma_i^{-2}x_i^2\right)d^nx.$$In the $n=2$ case, polar coordinates write the above as$$\int_0^{2\pi}d\theta\frac{1}{\sqrt{\pi}\sigma_1\sigma_2}(\sigma_1^{-2}+\sigma_2^{-2}+(\sigma_1^{-2}-\sigma_2^{-2})\cos2\theta)^{-3/2}.$$Define $a:=\frac{\sigma_1^{-2}-\sigma_2^{-2}}{\sigma_1^{-2}+\sigma_2^{-2}}$ so$$\int_0^{2\pi}(1+a\cos2\theta)^{-3/2}d\theta=4\int_0^{\pi/2}(1+a\cos2\theta)^{-3/2}d\theta=2\int_0^\pi(1+a\cos u)^{-3/2}du,$$which is already more than WA can handle.