3
$\begingroup$

Suppose you have a general $n$-th dimensional random Gaussian vector with probability distribution function $\mathcal{N}\left(\mathbf{x}|\boldsymbol{\mu},\boldsymbol{\Sigma}\right)$. What is the computational complexity (both space or time) of calculating the cumulative distribution function

$\int_{-\infty}^{y_1}dx_{1}\cdots\int_{-\infty}^{y_n}dx_{n}\mathcal{N}\left(\mathbf{x}|\boldsymbol{\mu},\boldsymbol{\Sigma}\right) $

?

This should be a function of $n$. If the suggested integration method is approximate (e.g., using monte-carlo methods), then the complexity may also depend on some required mean relative accuracy (std/mean) $\epsilon$.

Thanks in advance!

Improve this question
$\endgroup$
1

1 Answer 1

Reset to default
5
$\begingroup$

To some extent this depends on the model of computation. Some remarks:

  1. The problem must always involve an accuracy parameter $\epsilon$, since in general the answer will not be rational

  2. Even in the case of a standard 1-dimensional normal, the question is not so easy. I believe this can be done to accuracy $\epsilon$ in $\widetilde{O}(\log^2 (1/\epsilon))$ time -- i.e., essentially quadratic time in the number of digits of accuracy. The best reference here (I think) is Brent and Zimmerman's book "Modern Computer Arithmetic"; see their discussion of the erf and erfc functions.

  3. In the general $n$-dimensional case, I'm not sure whether you can get a $\mathrm{poly}(n, \log(1/\epsilon))$ running time. I'd possibly guess 'yes', but I don't know any reference. I would say that "in practice" the easiest solution would indeed be a naive Monte Carlo estimation. Assume first that you can generate standard Gaussians and do real arithmetic in constant time. Then you can get the answer (with high probability) to additive accuracy $\epsilon$ with $O(1/\epsilon^2)$ samples. Assuming you have the square-root of $\Sigma$, you could then probably get the answer in $O(n^2/\epsilon^2)$ time. (But then one technically has to worry about assumptions like generating standard Gaussians to sufficient accuracy, doing the real arithmetic, factorizing $\Sigma$ if necessary...)

Long story short, $\mathrm{poly}(n,1/\epsilon)$ time is surely okay, $\mathrm{poly}(n,\log(1/\epsilon))$ time I'm not sure about.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thank you for your detailed answer. However, I am more interested in relative accuracy $\epsilon$ rather than additive (this seems more relevant in high dimensions). This makes a large difference if we try to use monte-carlo. For concreteness, suppose we wish to calculate the positive orthant probability, where some components of $\boldsymbol{\mu} $ are strictly negative. In this case, since the relative area of the positive orthant in $\mathbb{R}^{n}$ is $2^{-n}$ , I guess we should get $O\left(c^{n}\epsilon^{-2}\right)$ complexity if we want relative accuracy (std/mean) of $\epsilon$. $\endgroup$
    – Daniel Soudry
    Commented Jan 3, 2014 at 16:33
  • 1
    $\begingroup$ Ah, good question. Okay, it occurs to me now that for relative accuracy $\epsilon$ you should be able to do it in time $\mathrm{poly}(n/\epsilon)$ using this paper of Cousins and Vempala (or those that it cites): arxiv.org/abs/1306.5829 It gives such an approximation for the Gaussian volume of any convex set, in fact. A small catch is that it requires some technical condition like the set containing the unit ball. However, from what I know of this subject, that is not a very strict requirement. It's more like a simple example of a possible requirement; really, the algorithm just $\endgroup$
    – Ryan O'Donnell
    Commented Jan 3, 2014 at 18:53
  • 1
    $\begingroup$ essentially needs to know a point inside the convex set to "get started". Since your set is an explicitly given intersection of n halfspaces, it should be straightforward to explicitly obtain such a point. In other words, I feel that it should be possible to work around the requirement and apply the Cousins-Vempala paper to your setting. $\endgroup$
    – Ryan O'Donnell
    Commented Jan 3, 2014 at 18:56
  • $\begingroup$ Thanks again - for this a very interesting paper, which indeed answers my general question. However, I do not think you can use this algorithm for my specific example. As far as I understand, the unit ball condition implies that the mean of the Gaussian is strictly contained within the set, so that the total probability mass in the set is non-vanishing with $n$. However, if the probability is exponentially vanishing (as in my example), then I think I can prove that the complexity is exponentially increasing in $n$ (as I wrote above) - no matter which sampling method is used. $\endgroup$
    – Daniel Soudry
    Commented Jan 4, 2014 at 13:04

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .