Let $X_i \sim \mathcal{N}(\mu_i, \sigma_i)$ be independent, normally-distributed random variables. Let
$$Y = a + b \max_i X_i$$
where $a \in \mathbb{R}$ and $b \in (0, 1)$. Which Gaussian distribution is closest to the distribution of $Y$? That is,
$$\operatorname*{argmin}_{\mu, \sigma} \text{KL}(F_Y, \mathcal{N}(\mu, \sigma))$$
where $F_Y$ and $f_Y$ denote the CDF and PDF of the distribution of $Y$, respectively. According to this answer we have
$$f_{\frac{Y-a}{b}}(y) = \left( \frac{1}{\sigma} \sum_i \frac{\phi \left( \frac{y-\mu_i}{\sigma} \right) }{\Phi \left( \frac{y-\mu_i}{\sigma} \right)} \right) \left( \prod_i \Phi \left( \frac{y-\mu_i}{\sigma} \right) \right)$$
where $\Phi$ and $\phi$ denote the CDF and PDF of the standard normal distribution, respectively. Unfortunately, the integral needed to compute the KL divergence (before taking its gradient to find the argmin) seems difficult. Does this problem have an analytic solution? If not, what's the best way to compute it numerically?