A process of repeated convolution and conditioning and the resulting sequence of probability distributions

Question

I am interested in the following procedure that yields a sequence $D_1,D_2,\ldots,$ of probability distributions over $\mathbb{R}^n$.

Let $D_1$ be the $n$-dimensional Gaussian distribution with covariance $I_n/n$ (so that the expected squared norm of a vector is $1$) conditioned on the norm of the result being at most $1$. So, $D_1$ is a Gaussian with its tail chopped off. (The starting distribution is probably not too important.)

Then, for $i \geq 1$, let $D_{i+1}$ be the distribution obtained by convolving $D_i$ with itself and then conditioning on the result having norm at most $1$.

Equivalently, if $f_i$ is the PDF of $D_i$, then $f_{i+1} = (f_i \circ f_i) \cdot 1_{\|x\| \leq 1}/a_i$, where $a_i = \Pr_{x,y \sim D_i}[\|x+y\| \leq 1]$ is chosen so that $f_{i+1}$ is a probability distribution.

---

I think this is a relatively natural process (and it arises in an important context, as I discuss briefly below), but I'm having quite a hard time actually understanding how to reason about it. And, my attempts to search the literature have failed because I have no idea what to search for :).

So, I think what I'm really just asking for is for any help at all in understanding such a process (or a reference to work that has considered such things). (Of course, there are many natural specific questions that one can ask about this process. E.g., does it converge to some distribution, and if so, what is the limiting distribution?)

For my application, I am actually mostly interested in the $a_i$. It is clear that for all $i$, $a_i \geq \Pr_{x,y \sim S^{n-1}}[\|x+y\| \leq 1] \approx (3/4)^{n/2}$, i.e., $a_i$ is always at least the probability that two random vectors on a sphere sum to a shorter vector. However, it seems that this bound might not actually be tight. E.g., $a_1 \approx (\sqrt{e}/2)^{n/2} \approx 0.82^{n/2}$. Perhaps all of the $a_i$ are exponentially larger than this lower bound of $(3/4)^{n/2}$?

If the bound $(3/4)^{n/2}$ is in fact exponentially loose (for all $i$---or even just for, say, $i \lesssim \log n$), this could imply faster algorithms for computational lattice problems, which would have major implications for the security of lattice-based cryptography(!). For this reason, I would be happy even with a heuristic analysis of this process.

However, again, any help at all in understanding this process would be wonderful :).

(One curious fact that makes me think that $(4/3)^{n/2}$ is not the right answer: this procedure maps the uniform distribution on the sphere (which is in some sense the worst case distribution) to a distribution that puts significantly more mass on short vectors. In particular, if we pretend that $D_i$ is the uniform distribution over a sphere, then $D_{i+1}$ satisfies $\Pr_{x \sim D_{i+1}}[\|x\| = r] \approx r^{n}(1-r^2/4) \cdot (4/3)^{n/2}$. Notice that the latter distribution actually puts more mass on short vectors than even the Gaussian.)

---

@spencer-peters (my student) has numerically computed the limiting distribution in the 1-d case. Here's what it looks like in red, compared to a Gaussian in blue:

The limiting value of $a_i$ is roughly $0.8$.

Answer 1

$\newcommand{\R}{\mathbb R}\newcommand{\la}{\lambda}\newcommand{\Ga}{\Gamma}$We have \begin{equation*} f_{j+1}(x)=(f_j*f_j)(x)\,1(\|x\|<1)/A_j \end{equation*} for $j=0,1,\dots$ and $x\in\R^n$, where $A_j:=\int_{\|x\|<1}dx\,(f_j*f_j)(x)$.

For each $j$ and each $x$, the value of $(f_j*f_j)(x)$ is an $n$-fold integral.

If $f_0$ is spherically symmetric, then so will be $f_j$ for each $j$.

Moreover, then, for $n\ge2$, the $n$-fold integral can be reduced to a double one.

Indeed, in what follows, do assume that $f_0$ is spherically symmetric. Assume also that $n\ge2$. Then for each $j=0,1,\dots$ there is a measurable function $g\colon[0,\infty)\to[0,\infty)$ such that \begin{equation*} f_j(x)=g_j(\|x\|) \end{equation*} for all $x\in\R^n$. So, \begin{equation*} (f_j*f_j)(x)=(Gg_j)(\|x\|), \end{equation*} where \begin{equation*} (Gg)(u):=\int_0^1\int_0^1\la_n(E_{u,a,b,da,db})g(a)g(b), \end{equation*} where $\la_n$ is the Lebesgue measure over $\R^n$, \begin{equation*} E_{u,a,b,da,db}:=\{y\in\R^n\colon a<\|y\|<a+da,\,b<\|ue_1-y\|<b+db\}, \end{equation*} and $e_1:=(1,0,\dots,0)\in\R^n$.

Shaded in the picture below is the small approximate parallelogram that is the intersection of the set $E_{u,a,b,da,db}$ (with $(u,a,b,da,db)=(6,5,4,0.4,0.4)$) with the "upper" half-plane of the 2D plane through the points $0$, $ue_1$, and $y$, where $y$ is such that $\|y\|=a$ and $\|ue_1-y\|=b$; the "upper" half-plane is bounded by the line through $0,ue_1$ and contains the point $y$. The value of $r$ is the distance from the point $y$ to the line through $0,ue_1$.

For infinitesimal $da$ and $db$, the area of the parallelogram is $dA=da\,db\,\sin t$, where $t$ is the angle between the vectors $y$ and $ue_1-y$. Also, expressing the area of the triangle $0(ue_1)y$ in two different ways, we see that $r=\frac{ab}u\,\sin t$. Also, by the cosine theorem, $\cos t=\frac{a^2+b^2-u^2}{2ab}$ and hence $\sin t=\frac1{2ab}\,\sqrt{4a^2b^2-(a^2+b^2-u^2)^2}$. We also have the triangle inequalities $|a-b|\le u\le a+b$. So, \begin{equation*} \begin{aligned} &\la_n(E_{u,a,b,da,db})\\ &=dA\,c_{n-2}r^{n-2}\,1(|a-b|\le u\le a+b) \\ &=\frac{c_{n-2}}{2^{n-1}}\frac{da\,db}{abu^{n-2}}\,(4a^2b^2-(a^2+b^2-u^2)^2)^{(n-1)/2} \,1(|a-b|\le u\le a+b), \end{aligned} \end{equation*} where $c_k=\dfrac{2\pi^{(k+1)/2}}{\Ga((k+1)/2)}$ is the "length/surface area" of the $k$-dimensional unit sphere $\mathbb S_k$ in $\R^{k+1}$, so that $c_0=2$, $c_1=2\pi$, etc. So, for a certain positive real constant $b_n$ depending only on $n$ and for all $u\in(0,1]$, \begin{equation*} (Gg)(u)=\frac{b_n}{u^{n-2}} \int_0^1\frac{da\,g(a)}a\, \\ \times\int_{\max(a,u-a)}^{\min(1,u+a)} \frac{db\,g(b)}b\,(4a^2b^2-(a^2+b^2-u^2)^2)^{(n-1)/2}. \end{equation*}

Finally, for $j=0,1,\dots$ and $u\in[0,1]$, \begin{equation*} g_{j+1}(u)=(Hg_j)(u), \end{equation*} where \begin{equation*} (Hg)(u):=(Gg)(u)\Big/A(g) \end{equation*} and \begin{equation*} A(g):=\int_{\|x\|<1}dx\,(Gg)(\|x\|)=c_{n-1}\int_0^1 du\,u^{n-1}(Gg)(u). \end{equation*}

The case $n=1$ is exceptional and easier.

Below are the graphs $\{(u,g_j(u))\colon0<u\le1\}$ for constant $g_0$, for $j=1$ (red), $j=2$ (green), $j=3$ (blue), and for $n\in\{1,2,3,5,10,20,100\}$ (the amount of calculations is about the same for all these $n$):

It appears that the $g_j$'s converge very fast (strangely enough, the fastest convergence seems to be for $n=5$, as compared with the other $n\in\{1,2,3,5,10,20,100\}$).