# 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 – user193726 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. – user193726 Jul 20 '18 at 14:56
• @UthsavChitra Unfortunately, not. – user193726 Jul 1 '19 at 9:06

You can use the law of large numbers to approximate your expectation pretty easily.

Edit: Analyticaly you can multiply a bunch of normal cdf evaluations together. For $$i > 0$$

\begin{align*} \mathbb{P}\left( {I = i} \right) &= \mathbb{P}\left( {{X_i} = \mathop {\max {X_j}}\limits_{j = 1,n} } \right) \\ &= \prod_{j \ne i} \mathbb{P}\left( {{X_i} > X_j } \right) \\ &= \prod_{j \ne i} \mathbb{P}\left( Y_j \le 0 \right) \end{align*}

where $$Y_j \sim \text{Normal}(\mu_j - \mu_i, \Sigma_j + \Sigma_i)$$.

If we code $$i=0$$ as the event that no data point is clearly a maximum, then $$\mathbb{P}\left( I = 0 \right) = 1 - \sum_{k=1}^n \mathbb{P}\left( I = k \right)$$

• 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 }}$ – user193726 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 – user193726 Jul 20 '18 at 19:28
• Thanks. Yes, my question arises from a specific problem where $n$ is given by $n = \prod\limits_{k = 1}^d {{n_k}}$ where the ${{n_k}}$'s are positive integers and $d$ is the problem dimension. In practice, $d$ can be up to, say, 30, so that $n$ can be extremely large (typically ${n_k} \sim 100$). That's why I'm looking for analytic approximation formulae rather than MC but for the time being that's the best we can do for small values of $d$ – user193726 Jul 23 '18 at 10:41
• @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