Approximating the mathematical expectation of the argmax of a Gaussian random vector

Let $X = \left( {{X_1},...,{X_n}} \right) \sim \mathcal{N}\left( {{\mathbf{\mu }},{\mathbf{\Sigma }}} \right)$ be a Gaussian random vector and $I = \mathop {\arg \max }\limits_{i = 1,n} {X_i}$.

$I$ has probability mass function

$\mathbb{P}\left( {I = i} \right) = \mathbb{P}\left( {{X_i} = \mathop {\max {X_j}}\limits_{j = 1,n} } \right) = \mathbb{P}\left( {{X_i} - \mathop {\max {X_j}}\limits_{j \ne i} > 0} \right)$

and mathematical expectation

$\mathbb{E}I = \sum\limits_{i = 1}^n {i\mathbb{P}\left( {I = i} \right)}$

Generally speaking, for large $n$ and arbitrary covariance matrix ${\mathbf{\Sigma }}$ , computing $\mathbb{E}I$ is very difficult because it requires the numerical evaluation of high-dimensional normal orthant integrals. So, apart from the IID and INID cases with a diagonal covariance matrix ${\mathbf{\Sigma }}$, banded covariance matrices and degenerate cases such as ${\mu _j} \gg {\mu _{i \ne j}}$ , under which conditions on ${\mathbf{\Sigma }}$ (e.g. correlation decay) can we get simple, easy-to-evaluate numerical approximations to $\mathbb{E}I$ (and $\mathbb{V}I$ as well)?

The covariance matrices ${\mathbf{\Sigma }}$ I'm interested in look like this:

Related questions:

• By "argmax" are you referring to what we call the "maximum order statistic" if we view the random vector as a sample, or to something else? – Alecos Papadopoulos Jul 20 '18 at 14:28
• @AlecosPapadopoulos Yes the argmax is the maximum order statistic also noted ${X_{(n)}}$ in extreme value theory – Fabrice Pautot Jul 20 '18 at 14:31
• It is not at all certain that we know of conditions that will allow for "easy numerical approximations" of the $n$-dimensional normal integral with correlation. – Alecos Papadopoulos Jul 20 '18 at 14:46
• @AlecosPapadopoulos ...Hence my question. For instance, I know from some papers in extreme value theory that if the correlations/covariances decay sufficiently fast, we can get simplified formulae for the maximum ${X_{(n)}}$. Therefore, we might get simplified formulae for the argmax as well, but I've not been able to find a single paper about this problem. Besides, sorry, ${X_{(n)}}$ is the maximum, not the argmax which is the index of the maximum. – Fabrice Pautot Jul 20 '18 at 14:56
• @UthsavChitra Unfortunately, not. – Fabrice Pautot Jul 1 '19 at 9:06

You can use the law of large numbers to approximate your expectation pretty easily. Here's some R code:

library(MASS)

getExpectationI <- function(numSamps, mu, Sigma){
samps <- mvrnorm(numSamps, mu, Sigma)
Is <- apply(samps, 1, which.max)
Ibar <- mean(Is)
SE.I <- sd(Is)
return(list(lower = Ibar-2*SE.I, est = Ibar, upper = Ibar + 2*SE.I))
}

numSamps <- 1000
mu <- as.matrix(c(1,2,3))
Sigma <- matrix(c(1,.5,.5,.5,1,.5,.5,.5,1),nrow=3)
getExpectationI(numSamps, mu, Sigma)


For the samples $\mathbf{X}^1,\mathbf{X}^2,\ldots,\mathbf{X}^{1000}$, this will calculate $$1000^{-1}\sum_{k=1}^{1000} I^k$$ where $I^k = \sum_{n=1}^3 n \mathbb{1}(\max_m X_m = X_n)$ is the index of the maximum of the $k$th $3$-d sample. For this example mean/covariance pair, it looks like $\mathbb{E}[I] \approx 2.802$.

• Yes, sampling from $X$ is the easy way. But what if $n=81$ (as above), $n=1000$ or $n=100000$ (my own $n$'s are virtually unbounded)? Can we efficiently sample (e.g. in polynomial time) an arbitrary large multivariate normal distribution (I should know, but I don't) and how many samples do we need in order to properly estimate $\mathbb{E}I$? If there is an efficient Monte Carlo method for arbitrary $n$ and ${\mathbf{\Sigma }}$ then the problem is solved. Otherwise, we'd like to get efficient, analytic approximations depending on the particular structure of ${\mathbf{\Sigma }}$ – Fabrice Pautot Jul 20 '18 at 19:09
• Both of those sound like separate questions. If you want to simulate quickly for higher dimensional data, you might have to write more bespoke code and introduce parallelization. Required sample size calculations can be done using the CLT. I can add that in to the answer shortly if you’re interested. – Taylor Jul 20 '18 at 19:17
• ...sure I am. As stated in the question (and other related questions such as stats.stackexchange.com/questions/326320/…, I'm considering the asymptotic case. Is there a polynomial time Monte Carlo algo for arbitary ${\mathbf{\Sigma }}$? If not, is there a polytime algo is some special cases? I'll add this in the question if you don't mind – Fabrice Pautot Jul 20 '18 at 19:28
• @FabricePautot I didn't realize you meant asymptotic when you said "large $n$". I also didn't realize this question was about coding up a sampler in parallel. Simulating iid normals in parallel is easy, as long as you don't screw up the seeds, but multiplying this by a scale matrix (and getting the scale matrix in the first place) might be difficult depending on $\Sigma$. I am not an expert in this area, so I can't make any broad guarantees on the computing time for arbitrary covariance matrices. Are you interested in a specific problem? – Taylor Jul 20 '18 at 20:28
• @FabricePautot "How many samples do we need in order to properly estimate $\mathbb{E}(I)$?" We could say that the sample distribution (for the mean of a sample from $I$) can be well approximated by a normal distribution. The standard deviation of the sample distribution will depend on the standard deviation of the distribution of $I$. This standard deviation will be in the worst case equal to $0.5 \cdot n$ (when the distribution of $I$ concentrates on the ends $P(I=1)=(I=n)=0.5$). You can use this to compute the desired necessary sample size for the Monte Carlo simulation. – Sextus Empiricus Sep 2 '19 at 11:38