Depth Degeneracy in Neural Networks: Vanishing Angles in Fully Connected
  ReLU Networks on Initialization
Despite remarkable performance on a variety of tasks, many properties of deep
neural networks are not yet theoretically understood. One such mystery is the
depth degeneracy phenomenon: the deeper you make your network, the closer your
network is to a constant function on initialization. In this paper, we examine
the evolution of the angle between two inputs to a ReLU neural network as a
function of the number of layers. By using combinatorial expansions, we find
precise formulas for how fast this angle goes to zero as depth increases. These
formulas capture microscopic fluctuations that are not visible in the popular
framework of infinite width limits, and leads to qualitatively different
predictions. We validate our theoretical results with Monte Carlo experiments
and show that our results accurately approximate finite network behaviour.
\review{We also empirically investigate how the depth degeneracy phenomenon can
negatively impact training of real networks.} The formulas are given in terms
of the mixed moments of correlated Gaussians passed through the ReLU function.
We also find a surprising combinatorial connection between these mixed moments
and the Bessel numbers that allows us to explicitly evaluate these moments.
Cameron Jakub, Mihai Nica
Depth Degeneracy in Neural Networks: Vanishing Angles in Fully Connected ReLU Networks on Initialization

###### Abstract

Stacking many layers to create truly deep neural networks is arguably what has led to the recent explosion of these methods. However, many properties of deep neural networks are not yet understood. One such mystery is the depth degeneracy phenomenon: the deeper you make your network, the closer your network is to a constant function on initialization. In this paper, we examine the evolution of the angle between two inputs to a ReLU neural network as a function of the number of layers. By using combinatorial expansions, we find precise formulas for how fast this angle goes to zero as depth increases. The formulas are given in terms of the mixed moments of correlated Gaussians passed through the ReLU function. We also find a surprising combinatorial connection between these mixed moments and the Bessel numbers.

**Keywords:** infinite limits of deep neural networks

## 1 Introduction

The idea of stacking many layers to make truly _deep_ neural networks (DNNs) is what arguably led to the neural net revolution in the 2010s. Indeed, from a function-space point of view, it is known that depth exponentially improves expressibility (Poole et al., 2016; Eldan and Shamir, 2015). However, an important but less well known fact is that under standard initialization schemes, deep neural networks become more and more _degenerate_ as depth gets larger and larger. One sense in which this happens is the phenomenon of vanishing and exploding gradients (Hanin, 2018). Another sense in which networks become degenerate is that a neural network gets closer and closer to a (random) constant function i.e. the network sends all inputs to the same output and cannot distinguish input points. This phenomenon seems to have been independently rediscovered from different points of view by several authors (Avelin and Karlsson, 2022; Dherin et al., 2022; Li et al., 2022; Hayou et al., 2019; Schoenholz et al., 2016).

One method already proposed to deal with the degeneracy phenomenon is the idea of activation function _shaping_. In particular, Li et al. (2022), showed that by rescaling the non-linear activation function (i.e. using Leaky ReLUs with leakiness depending on network depth) can lead to a non-trivial angle between inputs. This result essentially by rescaling the non-linear activation function used in the network in such a way that the Markov chain for
the evolution of the correlation between two inputs rescales to a non-trivial SDE. However, a detailed analysis of the evolution of the angle \(\theta\)_without_ any scaling (e.g. using an ordinary ReLU in all layers) remained an outstanding problem. This is the gap we fill in this article.

### Main Results for the Angle Process \(\theta^{\ell}\)

In this paper we examine the evolution of the _angle_\(\theta^{\ell}\) between two arbitrary inputs \(x_{\alpha},x_{\beta}\in\mathbb{R}^{n_{in}}\) after passing through \(\ell\) layers of a fully connected ReLU network (a.k.a. a multi-layer perceptron) on initialization. The angle is defined in the usual way by the inner product between two vectors in \(\mathbb{R}^{n_{\ell}}\).

\[\cos\left(\theta^{\ell}\right):=\frac{\langle F^{\ell}(x_{\alpha}),F^{\ell}(x _{\beta})\rangle}{\|F^{\ell}(x_{\alpha})\|\|F^{\ell}(x_{\beta})\|},\]

where \(n_{\ell}\) is the width (i.e. number of neurons) of the \(\ell\)-th layer and \(F^{\ell}:\mathbb{R}^{n_{in}}\to\mathbb{R}^{n_{\ell}}\) is the (random) neural network function mapping input to the post-activation logits in layer \(\ell\) on initialization. We assume here that the initialization is done with appropriately scaled independent Gaussian weights so that the network is on the "edge of chaos" (Hayou et al., 2019; Schoenholz et al., 2016), where the variance of each layer is order 1. (See Table 2 for our precise definition of the fully connected ReLU neural network).

With this setup, since the effect of each layer is independent of everything previous, \(\theta^{\ell}\) can be thought of as a Markov chain evolving as layer number \(\ell\) increases. As expected by the aforementioned "large depth degeneracy" phenomenon, we observe that the angle concentrates \(\theta^{\ell}\to 0\) as \(\ell\to\infty\) (see Figure 1 for an illustration). This indicates that the hidden layer representation of _any_ two inputs in a deep neural network become closer and closer to co-linear as depth increases.

In this paper, we obtain a simple, yet remarkably accurate, approximation for the evolution of \(\theta^{\ell}\) as a function of \(\ell\) that captures precisely how quickly this degeneracy happens for small angles \(\theta^{\ell}\) and large layer widths \(n_{\ell}\).

**Approximation 1**: _For small angles \(\theta^{\ell}\ll 1\) and large layer width \(n_{\ell}\gg 1\), the angle \(\theta^{\ell+1}\) at layer \(\ell+1\) is well approximated by_

\[\ln\sin^{2}(\theta^{\ell+1})\approx\ln\sin^{2}(\theta^{\ell})-\frac{2}{3\pi} \theta^{\ell}-\rho(n_{\ell}), \tag{1}\]

_where \(\rho(n_{\ell})\) is a constant which depends on the width \(n_{\ell}\) of layer \(\ell\), namely:_

\[\rho(n):=\ln\left(\frac{n+5}{n-1}\right)-\frac{10n}{\left(n+5\right)^{2}}+ \frac{6n}{\left(n-1\right)^{2}}=\frac{2}{n}+\mathcal{O}\left(n^{-2}\right). \tag{2}\]

Figure 1 illustrates how well this prediction matches Monte-Carlo simulations of \(\theta^{\ell}\) sampled from real networks. Also illustrated is the _infinite width_ prediction for \(\theta^{\ell}\) (discussed in Appendix A.9) which is less accurate at predicting finite width network behaviour than our formula, due to the \(n_{\ell}^{-1}\) effects that our formula captures in the term \(\rho(n_{\ell})\) but are not present in the infinite width formula. Approximation 1 is a simple corollary to the mathematically rigorous statement, Theorem 1, for the mean and variance of the random variable \(\ln\sin^{2}(\theta^{\ell})\).
#### 1.1.1 Consequences and Comparison to Previous Work

The simple Approximation 1 predicts that \(\theta^{\ell}\to 0\)_exponentially fast_ in \(\ell\) due to the term \(\rho(n)\); it predicts

\[\theta^{\ell}\leq\exp\left(-\sum_{i=1}^{\ell}\rho(n_{i})\right)=\exp\left(-\sum_ {i=1}^{\ell}\frac{2}{n_{i}}+\mathcal{O}(n_{i}^{-2})\right).\]

(Note that the exponential behaviour vanishes when \(n_{\ell}\to\infty\) with \(\ell\) fixed). In contrast to this prediction, an analysis using only expected values or equivalently working in the infinite-width \(n\to\infty\) limit predicts that \(\theta^{\ell}\to 0\) like \(\ell^{-1}\), which is qualitatively very different! The prediction of this rate was first demonstrated under the name "edge of chaos" (Hayou et al., 2019; Schoenholz et al., 2016) and again in Roberts et al. (2022); Hanin (2022). These earlier works working with the correlations showed that \(1-\cos(\theta_{\ell})\to 0\) like \(\ell^{-2}\), which is equivalent to \(\theta_{\ell}\to 0\) like \(\ell^{-1}\) by Taylor series expansion.

We can recover this infinite width prediction by replacing \(\rho(n)\) with \(0\) in our update rule (1). Exponentiating both sides and using \(\sin(\theta)\approx\theta,e^{\theta}\approx 1+\theta\) for \(\theta\ll 1\), Approximation 1 becomes \((\theta^{\ell+1})^{2}\approx\left(\theta^{\ell}\right)^{2}(1-\frac{2}{3\pi} \theta^{\ell})\), which is equivalent to the result of Proposition C.1 of

Figure 1: In the above figure, we feed 2 inputs with initial angle \(\theta^{0}=0.3\) and plot the value of \(\ln(\sin^{2}\theta^{\ell})\) at layer \(\ell\) of the network on initialization. The experiment is repeated 1000 times with independent initialization and a network width of \(n_{\ell}=256\); the 99% confidence interval for the Monte Carlo estimated mean \(\hat{\mathbf{E}}[\ln(\sin^{2}\theta^{\ell})]\) is displayed. On top of this, we plot our simple prediction method as in (1), as well as the infinite width prediction as in (12).
Hanin (2022). In that paper, the rule was derived directly from the infinite width update rule for \(\cos(\theta)\), 1 where they also show that this leads to \(\theta^{\ell}\approx\frac{1}{\ell}\) as \(\ell\to\infty\).

Footnote 1: Note that the recursion in Proposition C.1 is written in terms of the variable \(\epsilon_{\alpha\beta}^{\ell}\) instead of \(\theta\). The connection to our work is \(\epsilon_{\alpha\beta}^{\ell}=\frac{1-\cos(\theta^{\ell})}{2}\approx(\theta^{ \ell})^{2}/4\) as \(\theta\to 0\))

#### 1.1.2 More detailed results for the mean and variance

Approximation 1 comes from a simplification of more precise formulas for the mean and variance of the random variable \(\ln(\sin^{2}(\theta^{\ell}))\), which are stated in Theorem 1 below.

**Theorem 1** (Formula for mean and variance in terms of J functions): _Conditionally on the angle \(\theta^{\ell}\) in layer \(\ell\), the mean and variance of \(\ln\sin^{2}(\theta^{\ell+1})\) obey the following limit as the layer width \(n_{\ell}\to\infty\)_

\[\mathbf{E}[\ln\sin^{2}(\theta^{\ell+1})]= \mu(\theta^{\ell},n_{\ell})+\mathcal{O}(n_{\ell}^{-2}),\quad \mathbf{Var}[\ln\sin^{2}(\theta^{\ell+1})]=\sigma^{2}(\theta^{\ell},n_{\ell}) +\mathcal{O}(n_{\ell}^{-2}) \tag{3}\] \[\mu(\theta,n):= \ln\left(\frac{(n-1)(1-4J_{1,1}^{2})}{4J_{2,2}-1+n}\right)+ \frac{4(J_{2,2}+1)}{n\left(\frac{4J_{2,2}-1}{n}+1\right)^{2}}\] (4) \[-\frac{4\left(8J_{1,1}^{2}J_{2,2}-8J_{1,1}^{4}+4J_{1,1}^{2}-8J_{1,1}J_{3,1}+J_{2,2}+1\right)}{n\left(1-\frac{1}{n}\right)^{2}\left(1-4J_{1,1}^{ 2}\right)^{2}}\] \[\sigma^{2}(\theta,n):= \frac{8n(J_{2,2}+1)}{(4J_{2,2}-1+n)^{2}}+\frac{8n(8J_{1,1}^{2}J_{ 2,2}-8J_{1,1}^{4}+4J_{1,1}^{2}-8J_{1,1}J_{3,1}+J_{2,2}+1}{(n-1)^{2}(1-4J_{1,1}^ {2})^{2}}\] (5) \[-\frac{16n(2J_{1,1}^{2}-4J_{1,1}J_{3,1}+J_{2,2}+1)}{(4J_{2,2}-1+n )(n-1)(1-4J_{1,1}^{2})},\]

_where \(\mathbf{E},\mathbf{Var}\) denote the conditional mean and variance of quantities in layer \(\ell+1\) given the value of \(\theta^{\ell}\) in the previous layer and \(J_{a,b}:=J_{a,b}(\theta^{\ell})\) are the joint moments of correlated Gaussians passed through the ReLU function \(\varphi(x)=\max\{x,0\}\), namely_

\[J_{a,b}(\theta):=\mathbf{E}_{G,\hat{G}}[\varphi^{a}(G)\varphi^{b}(\hat{G})], \tag{6}\]

_where \(G,\hat{G}\) are marginally \(\mathcal{N}(0,1)\) random variables with correlation \(\mathbf{E}[G\hat{G}]=\cos(\theta)\)._

The joint moments \(J_{a,b}(\theta)\) are discussed in detail in Section 3. A new combinatorial method of computing these moments is presented, which is used to give an explicit formula is given for these joint-moments, which is presented in Theorem 4. Using the explicit formula for \(J_{a,b}\), the result of Theorem 1 can be used to obtain useful asymptotic formula for \(\mu\) and \(\sigma\), as in the following corollary.

**Corollary 2** (Small \(\theta\) asymptotics for mean and variance): _Conditionally on the angle \(\theta^{\ell}\) in layer \(\ell\), the mean and variance of \(\ln\sin^{2}(\theta^{\ell+1})\) obey the following limit as the layer
_width \(n_{\ell}\to\infty\)_

\[\mathbf{E}[\ln\sin^{2}(\theta^{\ell+1})]= \mu(\theta^{\ell},n_{\ell})+\mathcal{O}(n_{\ell}^{-2}),\quad \mathbf{Var}[\ln\sin^{2}(\theta^{\ell+1})]=\sigma^{2}(\theta^{\ell},n_{\ell})+ \mathcal{O}(n_{\ell}^{-2}) \tag{7}\] \[\mu(\theta,n)= \ln\sin^{2}\theta-\frac{2}{3\pi}\theta-\rho(n)-\frac{8\theta}{15 \pi n}-\left(\frac{2}{9\pi^{2}}-\frac{68}{45\pi^{2}n}\right)\theta^{2}+ \mathcal{O}(\theta^{3})\] (8) \[\sigma^{2}(\theta,n)= \frac{8}{n}-\frac{64}{15\pi}\frac{\theta}{n}-\left(8+\frac{296}{4 5\pi}\right)\frac{\theta^{2}}{n}+\mathcal{O}\left(\theta^{3}\right), \tag{9}\]

_where \(\rho(n)\) is as defined in (2)._

To derive Approximation 1 from Theorem 1, we simply keep only the first few terms of the series expansion (8), and then also completely drop the variability, essentially approximating \(\sigma^{2}(\theta^{\ell},n)\approx 0\) (Note that in reality \(\sigma^{2}(\theta,n)\approx 8/n\) from (9)). Therefore Approximation 1 is a greatly simplified consequence of Theorem 1.

Moreover, our derivation shows that \(\ln\sin^{2}(\theta^{\ell})\) can be expressed in terms of averages over \(n\) pairs of independent Gaussian variables (see (14-16)). Thus, by central-limit-theorem type arguments, one would expect the following approximation by Gaussian laws which also accounts for the variability of \(\ln\sin^{2}(\theta)\) using our calculated value for the variance.

**Approximation 2**: _Conditional on the value of \(\theta^{\ell}\), the angle at layer \(\ell+1\) is well approximated by a Gaussian random variable_

\[\ln\sin^{2}(\theta^{\ell+1})\stackrel{{ d}}{{\approx}}\mathcal{N }(\mu(\theta^{\ell},n_{\ell}),\sigma^{2}(\theta^{\ell},n_{\ell})), \tag{10}\]

_where \(\mu,\sigma^{2}\) are as in Theorem 2._

Figure 2: Plots comparing the functions \(\mu(\theta,n)\) and \(\sigma^{2}(\theta,n)\) to simulated neural networks. The linear approximation of \(\mu\), used to create Approximation 1 is also displayed. Confidence bands are constructed by randomly initializing 10000 neural networks with layer width \(n_{\ell}=1024\), and a range of 100 initial angles \(0.005\leq\theta^{\ell}\leq 0.8\). We study \(\theta^{\ell+1}\) and use the simulations to construct 99% confidence intervals for a) \(\mathbf{E}\left[\ln(\sin^{2}(\theta^{\ell}))-\ln(\sin^{2}(\theta^{\ell+1}))\right]\) and b) \(\mathbf{Var}\left[\ln(\sin^{2}(\theta^{\ell+1}))\right]\).
We find that the normal approximation (10) matches simulated finite neural networks remarkably well; see Monte-Carlo simulations from real networks in Figure 3.

The methods we use to obtain these approximations are quite flexible. For example, more accurate approximations can be obtained by incorporating higher moments \(J_{a,b}(\theta)\) (see Section 2 for a discussion). We also believe that it should be possible to extend these methods to other non-linearities beyond ReLU and more complicated neural network architectures through the same basic principles we introduce here.

### J functions and infinite width limits

In 2009, Cho and Saul (Cho and Saul (2009)) introduced the \(p\)-th moment for correlated ReLU-Gaussians, which they denoted with the letter \(J\),

\[J_{p}(\theta):=2\pi\mathbf{E}\left[\varphi^{p}(G)\varphi^{p}(\hat{G})\right], \tag{11}\]

where \(p\in\mathbb{N}\), and \(\varphi(x)=\max\{x,0\}\) is the ReLU function, \(G,\hat{G}\in\mathbb{R}\) are marginally two standard \(\mathcal{N}(0,1)\) Gaussian random variables with correlation \(\mathbf{Cov}(G,\hat{G})=\cos(\theta)\). This quantity has found numerous applications for infinite width networks. One simple application of \(J_{1}\) appears when in the infinite width approximation for \(\cos(\theta^{\ell})\), where \(\ell\) is fixed an we take the limit \(n_{1},n_{2},\ldots n_{\ell}\rightarrow\infty\) (see Appendix A.9 for a detailed derivation):

**Approximation 3**: _The infinite-width approximation for the angle \(\theta^{\ell+1}\) given \(\theta^{\ell}\) is_

\[\cos\left(\theta^{\ell+1}\right)=\frac{J_{1}(\theta^{\ell})}{\pi}=\frac{\sin( \theta^{\ell})+(\pi-\theta^{\ell})\cos(\theta^{\ell})}{\pi}. \tag{12}\]

Figure 3: Histograms of \(\ln(\sin^{2}(\theta^{\ell+1}))\) using 5000 simulated neural networks with an initial angle of \(\theta^{\ell}=0.1\), and layer width \(n_{\ell}=256\) for all layers \(1\leq\ell\leq 30\). We plot the simulated histograms as well as our predicted probability density function using Approximation 2 from (10) at layer 1 (left) and layer 30 (right). The predicted and empirical distribution are statistically indistinguishable according to a A Kolmogorov-Smirnov test, with \(p\) values \(0.913>0.05\) (left) and \(0.450>0.05\) (right).
The formula for \(J_{1}\) is the \(p=1\) case of a remarkable explicit formula for \(J_{p}\) derived by Cho and Saul (2009) namely,

\[J_{p}(\theta)=(-1)^{p}(\sin\theta)^{2p+1}\left(\frac{1}{\sin\theta}\frac{ \partial}{\partial\theta}\right)^{p}\left(\frac{\pi-\theta}{\sin\theta}\right).\]

This allows one to derive asymptotics of \(\theta^{\ell}\) in the infinite width limit. However, there are several limitations to this approach. Most important is that infinite width limit is not a good approximation when the network depth \(\ell\) is comparable to the network width \(n\)(Li et al. (2022)). The infinite width limit uses the law of large numbers to obtain (12), thereby discarding random fluctuations. For very deep networks, microscopic fluctuations (on the order of \(\mathcal{O}(1/n_{\ell})\)) from layer to layer can accumulate over \(\ell\) layers to give macroscopic effects. This is why the infinite width predictions for \(\theta^{\ell}\) are not a good match to the simulations in Figure 1; very deep networks are far from the infinite width limit in this case. See Figure 1 where the infinite width predicitons are compared to finite networks.

Instead, to analyze the evolution of the angle \(\theta^{\ell}\) more accurately, we need to do something more precise than the law of large numbers to capture the effect of these microscopic fluctuations. This is the approach we carry out in this paper. While the mean only depends on the \(p\)-th moment functions \(J_{p}\) from (11), these fluctuations depend on the _mixed_ moments, which we denote by \(J_{a,b}\) for \(a,b\in\mathbb{N}\) as follows2

Footnote 2: Note that compared to Cho and Saul’s definition for \(J_{p}\), we omit the factor of \(2\pi\) in our definition of \(J_{a,b}\). The factor of \(2\pi\) seems natural when \(a+b\) is even (like the case \(a=b=p\) that Cho-Saul considered), but when \(a+b\) is odd a different factor of \(2\sqrt{2\pi}\) appears! Therefore the factor of \(2\pi\) would confuse things in the general case (see Table 1). The correct translation between Cho-Saul \(J_{p}\) and our \(J_{a,b}\) is \(J_{p}=2\pi J_{p,p}\)

\[J_{a,b}(\theta):=\mathbf{E}\left[\varphi^{a}(G)\varphi^{b}(\hat{G})\right], \tag{13}\]

with \(G,\hat{G}\) again as in (11) are marginally \(\mathcal{N}(0,1)\) with correlation \(\cos(\theta)\). In Section 2.1 we carry out a detailed asymptotic analysis to write the evolution of \(\theta^{\ell}\) in terms of the mixed moments \(J_{a,b}\). In order to make useful predictions, one must also calculate a formula for \(J_{a,b}(\theta)\). Unfortunately, the method that Cho-Saul originally proposed for this does _not_ seem to work when \(a\neq b\). This is because that method used contour integrals, and relied on using certain trig identities which do not hold when \(a\neq b\). Instead, in Section 3, we introduce a new method, based on Gaussian integration by parts, to compute \(J_{a,b}\) for general \(a,b\) via a recurrance relation. By serendipity3, we find a remarkable combinatorial connection between \(J_{a,b}\) and the Bessel numbers (Cheon et al. (2013)), which allows one to find an explicit (albeit complicated) formula for \(J_{a,b}\) in terms of binomial coefficients. The formula for the first few functions are shown in Table 1, and the general explicit formula is presented in Theorem 4.

Footnote 3: This connection was first noticed by calculating the first few \(J\) functions, and then using the On-Line Encyclopedia of Integer Sequence to discover the connection to Bessel number [https://oeis.org/A001498](https://oeis.org/A001498)

### Outline

The two main contribution of this paper is to prove Theorem 1 for the evolution of \(\theta^{\ell}\) in terms of the mixed moments \(J_{a,b}\), and then seperat
for any of the mixed moments \(J_{a,b}\). Combined these allow for the explicit formula Corollary 2 and the useful simpler Approximations 1 and 2. See Figure 2 and Figure 3 for comparisons of these predictions to Monte-Carlo simulations. We also believe the methods proposed here are flexible enough to be modified to apply to non-linearities other than ReLU and to different neural network architectures beyond fully-connected networks in future work.

Section 2 contains the analysis of the angle process and predicted distribution of \(\ln(\sin^{2}(\theta^{\ell}))\) in deep ReLU networks. Section 2.1 covers our approximation of \(\mathbf{E}[\ln(\sin^{2}(\theta^{\ell+1}))]\) which leads to the rule for \(\theta^{\ell}\) as in equation (1), while Section 2.2 outlines our approximation for \(\mathbf{Var}[\ln(\sin^{2}(\theta^{\ell+1}))]\).

In Section 3 we cover the derivation of the explicit formula for the \(J\) functions. We state the main results of this section in Section 3.1, and cover the mathematical tools needed to solve the expectations using Gaussian integration by parts in Section 3.2. We develop the formula for \(J_{a,b}\) by first finding a recursive formula in Section 3.3, which reveals a connection between the \(J\) functions and the Bessel numbers. This recursion is studied to develop an explicit formula for \(J_{a,b}\) in Section 3.4.

## 2 ReLU Neural Networks on Initialization

In this section, we analyze ReLU neural networks and show how the the mixed moments \(J_{a,b}\) appear in evolution of the angle \(\theta^{\ell}\) on initialization. We define the notation we used for a fully connected ReLU neural network, along with other notations we will use in Table 2. Note that the factor of \(\sqrt{2/n_{\ell}}\) in our definition is implementing the so called He initialization (He et al., 2015), which ensures that \(\mathbf{E}[\|z^{\ell}\|^{2}]=\|x\|^{2}\) for all layers \(\ell\). This initialization is known to be the "critical" initialization for taking large limits of the network (Roberts et al., 2022; Hayou et al., 2019). Given this neural network, we wish to study the evolution

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline \(a\)\(b\) & \(0\) & \(1\) & \(2\) & \(3\) \\ \hline \(0\) & \(\pi-\theta\) & \(\cos\theta\)+\(1\) & \((\pi-\theta)\)+\(\sin\theta\cos\theta\) & \(2(\cos\theta\)+\(1)\)+\(\sin^{2}\theta\cos\theta\) \\ \hline \(1\) & & \(\sin\theta\)+\((\pi-\theta)\cos\theta\) & \((\cos\theta\)+\(1)^{2}\) & \(3(\pi-\theta)\cos\theta\)+\(\sin\theta\cos^{2}\theta\)+\(2\sin\theta\) \\ \hline \(2\) & & & \((\pi-\theta)(2\cos^{2}\theta\)+\(1)\)+\(3\sin\theta\cos\theta\) & \(3\cos\theta(\cos\theta\)+\(1)^{2}\)+\(2(\cos\theta\)+\(1)\)+\(\sin^{2}\theta\cos\theta\) \\ \hline \(3\) & & & & \((\pi-\theta)(6\cos^{2}\theta\)+\(9)\cos\theta\) \\  & & & & \(+5\sin\theta\cos^{2}\theta\)+\((6\cos^{2}\theta\)+\(4)\sin\theta\) \\ \hline \end{tabular}
\end{table}
Table 1: Table of formulas for the first few \(J\) functions. The normalizing constant appearing in all entries, either \(c_{0}=2\pi\), \(c_{1}=2\sqrt{2\pi}\) depending on the parity of \(a+b\), has been omitted, i.e. the table shows the value of \(c_{(a+b\bmod 2)}J_{a,b}(\theta)\). (Note that when \(a=0\), the appropriate convention of \(0^{0}=0\) is needed, see (25) for details). These generalize \(J_{p}(\theta)\) of (11) which appear on the diagonal of this table. Note that \(J_{a,b}=J_{b,a}\) so only upper triangular entries are shown. An explicit formula for all \(J_{a,b}\) is derived in Section 3.4.
of 2 inputs \(x_{\alpha}\) and \(x_{\beta}\) as they traverse through the layers of the network. Specifically, we wish to study how the angle \(\theta\) between the inputs changes as the inputs move from layer to layer.

The starting point for our calculation is to notice that because the weights are Gaussian, the values of \(\varphi_{\alpha}^{\ell+1},\varphi_{\beta}^{\ell+1}\) are jointly Gaussian given the vectors of \(\varphi_{\alpha}^{\ell},\varphi_{\beta}^{\ell}\). In fact, it turns out that by properties of Gaussian random variables, one only needs to know the values of the scalars \(\|\varphi_{\alpha}^{\ell}\|,\|\varphi_{\beta}^{\ell}\|\) and \(\theta^{\ell}\) to understand the full distribution of \(\varphi_{\alpha}^{\ell+1},\varphi_{\beta}^{\ell+1}\). (see Appendix A.5 for details) By using the positive homogeneity of the ReLU function \(\varphi(\lambda x)=\lambda\varphi(x)\) for \(\lambda>0\), we can factor out the effect of the norm of each vector in layer \(\ell\). After some manipulations, these ideas lead us to the following identities that are the heart of our calculations; a full derivation of these quantities are provided in Appendix A.5 and A.6.

\[\|\varphi_{\alpha}^{\ell+1}\|^{2} =\frac{\|\varphi_{\alpha}^{\ell}\|^{2}}{n_{\ell}}\sum_{i=1}^{n_{ \ell}}2\varphi^{2}(G_{i}), \|\varphi_{\beta}^{\ell+1}\|^{2} =\frac{\|\varphi_{\beta}^{\ell}\|^{2}}{n_{\ell}}\sum_{i=1}^{n_{\ell}}2 \varphi^{2}(\hat{G}_{i}), \tag{14}\] \[\langle\varphi_{\alpha}^{\ell+1},\varphi_{\beta}^{\ell+1}\rangle =\frac{\|\varphi_{\alpha}^{\ell}\|\|\varphi_{\beta}^{\ell}\|}{n _{\ell}}\sum_{i=1}^{n_{\ell}}2\varphi(G_{i})\varphi(\hat{G}_{i}),\] (15) \[\frac{\|\varphi_{\alpha}^{\ell+1}\|^{2}\|\varphi_{\beta}^{\ell+1} \|^{2}}{\|\varphi_{\alpha}^{\ell}\|^{2}\|\varphi_{\beta}^{\ell}\|^{2}}\sin^{2} (\theta^{\ell+1}) =\frac{2}{n_{\ell}^{2}}\sum_{i,j=1}^{n_{\ell}}\left(\varphi(G_{i}) \varphi(\hat{G}_{j})-\varphi(G_{j})\varphi(\hat{G}_{i})\right)^{2}, \tag{16}\]

\begin{table}
\begin{tabular}{|c|c|} \hline
**Symbol** & **Definition** \\ \hline \hline \(x\in\mathbb{R}^{n_{in}}\) & Input (e.g. training example) in the input dimension \(n_{in}\in\mathbb{N}\) \\ \hline \(\ell\in\mathbb{N}\) & Layer number. \(\ell=0\) is the input \\ \hline \(n_{\ell}\in\mathbb{N}\) & Width of hidden layer \(\ell\) (i.e. number of neurons in layer \(\ell\)) \\ \hline \(W^{\ell}\in\mathbb{R}^{n_{\ell+1}\times n_{\ell}}\) & Weight matrix for layer \(\ell\). Initialized with iid standard Gaussian entries \\  & \(W^{\ell}_{a,b}\sim\mathcal{N}(0,1)\) \\ \hline \(\varphi:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}\) & Entrywise ReLU activation function \(\varphi(x)_{i}=\varphi(x_{i})=\max\{x_{i},0\}\) \\ \hline \(z^{\ell}(x)\in\mathbb{R}^{n_{\ell}}\) & Pre-activation vector in the \(\ell^{\text{th}}\) layer for input \(x\) (a.k.a logits of layer \(\ell\)) \\  & \(z^{1}(x):=W^{1}x,\qquad z^{\ell+1}(x):=\sqrt{\frac{2}{n_{\ell}}}W^{\ell+1} \varphi(z^{\ell}(x))\). \\ \hline \(\varphi_{\alpha}^{\ell},\varphi_{\beta}^{\ell}\in\mathbb{R}^{n_{\ell}}\) & Post-activation vector on inputs \(x_{\alpha},x_{\beta}\) respectively \\  & \(\varphi_{\alpha}^{\ell}:=\varphi(z^{\ell}(x_{\alpha})),\qquad\varphi_{\beta}^{ \ell}:=\varphi(z^{\ell}(x_{\beta}))\) \\ \hline \(\theta^{\ell}\in[0,\pi]\) & Angle between \(\varphi_{\alpha}^{\ell}\) and \(\varphi_{\beta}^{\ell}\) defined by \(\cos(\theta^{\ell}):=\frac{\langle\varphi_{\alpha}^{\ell},\varphi_{\beta}^{ \ell}\rangle}{\|\varphi_{\alpha}^{\ell}\|\|\varphi_{\beta}^{\ell}\|}\) \\ \hline \(R^{\ell+1}\in\mathbb{R}\) & Shorthand for the ratio \(R^{\ell+1}:=\frac{\|\varphi_{\alpha}^{\ell+1}\|^{2}\|\varphi_{\beta}^{\ell+1} \|^{2}}{\|\varphi_{\alpha}^{\ell}\|^{2}\|\varphi_{\beta}^{\ell}\|^{2}}\) \\ \hline \end{tabular}
\end{table}
Table 2: Definition and notation used for fully connected ReLU neural networks
where \(G_{i}\), \(\hat{G}_{i}\) are all marginally \(\mathcal{N}(0,1)\), with correlation \(\mathbf{Cov}(G_{i},\hat{G}_{i})=\cos(\theta^{\ell})\) and independent for different indices \(i\). The identity in (16) is derived using the _determinant of the Gram matrix_ for vectors \(\varphi_{\alpha}^{\ell+1},\ \varphi_{\beta}^{\ell+1}\) (full derivation given in Appendix A.6). For convenience, we define the notation \(R^{\ell+1}=\frac{\|\varphi_{\alpha}^{\ell+1}\|^{2}\|\varphi_{\beta}^{\ell+1} \|^{2}}{\|\varphi_{\alpha}^{\ell}\|^{2}\|\varphi_{\beta}^{\ell}\|^{2}}\). Combining the equations (14) gives us a useful identity for \(R^{\ell+1}\):

\[R^{\ell+1}=\frac{4}{n_{\ell}^{2}}\sum_{i,j=1}^{n_{\ell}}\varphi^{2}(G_{i}) \varphi^{2}(\hat{G}_{j}). \tag{17}\]

Given some \(\theta^{\ell}\), we wish to predict the behaviour of \(\theta^{\ell+1}\). Rather than studying \(\theta^{\ell+1}\) directly, we instead study the quantity \(\ln(\sin^{2}(\theta^{\ell+1}))\). This allows us to use convenient approximations and identities for quantities we are interested in. (And indeed, a post-hoc analysis shows that as \(\theta\to 0\), the random variable \(\ln\sin^{2}(\theta^{\ell+1})\) has a _non-zero constant_ variance which depends only on \(n_{\ell}\). This is in contrast to \(\theta^{\ell}\) itself which has variance tending to _zero_. This is one reason why the Gaussian approximation for \(\ln\sin^{2}(\theta^{\ell})\in(-\infty,\infty)\) works well, whereas Gaussian approximations for \(\theta^{\ell}\) or \(\cos(\theta^{\ell})\in[-1,1]\) are less accurate.) We first derive a formula for \(\mathbf{E}\left[\ln(\sin^{2}(\theta^{\ell+1}))\right]\).

### Expected Value

In this section we show how to compute the expected value of \(\ln(\sin^{2}(\theta^{\ell}))\) in terms of the \(J\) functions as in Theorem 1. Firstly, we rewrite this expectation as the difference

\[\mathbf{E}\left[\ln(\sin^{2}(\theta^{\ell+1}))\right]=\mathbf{E}\left[\ln \left(R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right)\right]-\mathbf{E}\left[\ln \left(R^{\ell+1}\right)\right]. \tag{18}\]

The two random variables \(R^{\ell+1}\) and \(R^{\ell+1}\sin^{2}(\theta^{\ell+1})\) in (18) both have interpretations in terms of sums of Gaussians as in (17) and (16) which makes it possible to calculate their moments in terms of the \(J\) functions. To enable our use of the moments here, we use the following approximation of \(\ln(X)\) for a random variables \(X\), which is based on the Taylor expansion for \(\ln(1+x)=x-\frac{1}{2}x^{2}+\ldots\) (a full derivation is given in Appendix A.1):

\[\ln(X)=\ln(\mathbf{E}[X])+\frac{X-\mathbf{E}[X]}{\mathbf{E}[X]}-\frac{(X- \mathbf{E}[X])^{2}}{2\mathbf{E}[X]^{2}}+\epsilon_{2}\left(\frac{X-\mathbf{E}[ X]}{\mathbf{E}[X]}\right), \tag{19}\]

where \(\epsilon_{2}(x)\) is the Taylor remainder term in \(\ln(1+x)=x-\frac{x^{2}}{2}+\epsilon_{2}(x)\) and satisfies \(\epsilon_{2}(x)=\mathcal{O}(x^{3})\). Applying this approximation to the terms appearing on the right hand side of (18), and taking expected value of both sides, we obtain the estimates

\[\mathbf{E}\left[\ln\left(R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right)\right] =\ln\left(\mathbf{E}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1}) \right]\right)-\frac{\mathbf{Var}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1}) \right]}{2\mathbf{E}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right]^{2}}+ \mathcal{O}(n_{\ell}^{-2}),\] \[\mathbf{E}\left[\ln\left(R^{\ell+1}\right)\right] =\ln\left(\mathbf{E}\left[R^{\ell+1}\right]\right)-\frac{\mathbf{ Var}\left[R^{\ell+1}\right]}{2\mathbf{E}\left[R^{\ell+1}\right]^{2}}+ \mathcal{O}(n_{\ell}^{-2}).\]

To control the error here, we have used here the fact that \(R^{\ell+1}\) and \(R^{\ell+1}\sin^{2}(\theta^{\ell+1})\) can be written as averages over random variables as in (14 - 16). This allows us to show the 3rd
central moments for \(R^{\ell+1}\) and \(R^{\ell+1}\sin^{2}(\theta^{\ell+1})\) are \(\mathcal{O}(n_{\ell}^{-2})\); see Appendix A.1 for details. This approximation is convenient because we are able to calculate the values on the right hand side of the equations in terms of the moments \(J_{a,b}\) by expanding/taking expectations of the representations (14 - 16). The key quantities we calculate are

\[\mathbf{E}\left[R^{\ell+1}\right] =\frac{4J_{2,2}-1}{n_{\ell}}+1 \tag{20}\] \[\mathbf{E}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right] =\frac{(n_{\ell}-1)(1-4J_{1,1}^{2})}{n_{\ell}}\] (21) \[\mathbf{Var}\left[R^{\ell+1}\right] =\frac{4}{n_{\ell}}(J_{2,2}+1)+\frac{16}{n_{\ell}^{2}}\left(2J_{ 4,2}-\frac{5}{2}J_{2,2}+J_{2,2}^{2}+\frac{5}{8}\right)+\mathcal{O}\left(n_{ \ell}^{-3}\right)\] (22) \[\mathbf{Var}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right] =\frac{8\left(-8J_{1,1}^{4}+8J_{1,1}^{2}J_{2,2}+4J_{1,1}^{2}-8J_{ 1,1}J_{3,1}+J_{2,2}+1\right)}{n_{\ell}}+\mathcal{O}\left(n_{\ell}^{-2}\right), \tag{23}\]

where \(J_{a,b}=J_{a,b}(\theta^{\ell})\). These formulas are calculated in Appendix A.8 by a combinatorial expansion using the representations from (14-16). Combining these gives the result for \(\mu(\theta,n)\) in Theorem 1. Note that to obtain a more accurate approximation, we would simply include more terms in the variance expressions in (22, 23).

### Variance of \(\ln(\sin^{2}(\theta^{\ell+1}))\)

In this section we show how to compute the variance of \(\ln(\sin^{2}(\theta^{\ell}))\) in terms of the \(J\) functions as in Theorem 1. We can rewrite \(\mathbf{Var}[\ln(\sin^{2}(\theta^{\ell+1}))]\) in the following way:

\[\mathbf{Var}[\ln(\sin^{2}(\theta^{\ell+1}))]=\mathbf{Var}\left[ \ln\left(R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right)-\ln\left(R^{\ell+1} \right)\right] \tag{24}\] \[=\mathbf{Var}\left[\ln\left(R^{\ell+1}\sin^{2}(\theta^{\ell+1}) \right)\right]+\mathbf{Var}\left[\ln\left(R^{\ell+1}\right)\right]-2\mathbf{ Cov}\left(\ln\left(R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right),\ln\left(R^{ \ell+1}\right)\right).\]

We have now expressed this in terms of \(R^{\ell+1}\) and \(R^{\ell+1}\sin^{2}(\theta^{\ell+1})\) which will allow us to use identities as in (14 - 16) in our calculations. Appendix A.2 and Appendix A.3 cover the method used to approximate the unknown variance and covariance terms above. Once again, we control the error term arising from moments in the error term of the Taylor series by using representation as sums (14 - 16). We have already calculated most of the quantities on the right hand side already in our calculation for \(\mu(\theta,n)\), the only new term is

\[\mathbf{Cov}\left(R^{\ell+1}\sin^{2}(\theta^{\ell+1}),R^{\ell+1}\right)=\frac {1}{n_{\ell}}\left(16J_{1,1}^{2}-32J_{1,1}J_{3,1}+8J_{2,2}+8\right)+\mathcal{ O}\left(n_{\ell}^{-2}\right).\]

This is again computed by a combinatorial expansion of the sums (14-16). (Full calculation given in Appendix A.8). We now have solved for all of the functions needed to perform our approximation of \(\mathbf{Var}[\ln(\sin^{2}(\theta^{\ell+1}))]\). Putting it together, we end up with the expression for \(\sigma^{2}(\theta,n)\) as in (5). We compare the predicted probability distribution of \(\ln(\sin^{2}(\theta))\) using our formulas \(\mu(\theta,n)\) and \(\sigma^{2}(\theta,n)\) to empirical probability distributions in Figure 3.
## 3 Explicit formula for the mixed-moment J Functions

In this section we develop a combinatorial method that allows us compute exact formulas for the \(J\) functions. The method is to use Gaussian integration by parts to find a recurrence relationship between the moments \(J_{a,b}\), and then solve it explicitly. We begin by generalizing the definition of \(J_{a,b}\) from (13) to include \(a=0\) and/or \(b=0\) as follows. Let \(G\), \(W\) be _independent_\(\mathcal{N}(0,1)\) variables. Then, we define the functions \(J_{a,b}(\theta)\) as

\[J_{a,b}(\theta)=\mathbf{E}[G^{a}(G\cos\theta+W\sin\theta)^{b}\;1\{G>0\}\;1\{G \cos\theta+W\sin\theta>0\}], \tag{25}\]

where \(a,b\in\mathbb{N}\cup\{0\}\). Note that \(G\cos\theta+W\sin\theta=\hat{G}\) is marginally \(\mathcal{N}(0,1)\) and has correlation \(\cos(\theta)\) with \(G\), matching the original definition. The ReLU function satisfies the identity \(\varphi(x)^{a}=x^{a}1\{x>0\}\) for \(a\geq 1\), so (25) generalizes (13) to the case \(a=0\). We also note that \(J_{a,b}(\theta)=J_{b,a}(\theta)\) for all \(a,b\in\mathbb{N}\cup\{0\}\).

### Statement of Main Results and Outline of Method

By using the method of Gaussian integration by parts, we are able to derive recurrence relations for the \(J_{a,b}\). Since the definition of \(J_{a,b}\) involves the indicator function \(1\{G>0\}\), we must make sense of what the derivative of this function means for the purposes of integration by parts; see Section 3.2 where this is carried out. Then, by use of the generalized Gaussian integration by parts formula, we obtain the following recurrence relations for \(J_{a,b}\).

**Proposition 1** (Recurrence relations for \(J_{a,b}\)): _For \(a\geq 2\), the sequence \(J_{a,0}\) satisfies the recurrence relation:_

\[J_{a,0}(\theta)=(a-1)J_{a-2,0}(\theta)+\frac{\sin^{a-1}\theta\cos\theta}{c_{a \bmod 2}}(a-2)!!\,, \tag{26}\]

_where \(c_{0}=2\pi\), \(c_{1}=2\sqrt{2\pi}\). For \(a\geq 2\), and \(b\geq 1\), the collection, \(J_{a,b}\) satisfies the following two-index recurrence relation:_

\[J_{a,b}(\theta)=(a-1)J_{a-2,b}(\theta)+b\cos\theta J_{a-1,b-1}(\theta). \tag{27}\]

The same integration by parts technique that yields the recurrence relation also makes it easy to evaluate the first few \(J\) functions. They are as follows

**Proposition 2** (Explicit Formula for \(J_{0,0},J_{1,0},J_{1,1}\)): \(J_{0,0}\)_, \(J_{1,0}\), and \(J_{1,1}\) are given by_

\[J_{0,0}(\theta)=\frac{\pi-\theta}{2\pi},\qquad J_{1,0}(\theta)=\frac{1+\cos \theta}{2\sqrt{2\pi}},\qquad J_{1,1}(\theta)=\frac{\sin\theta+(\pi-\theta)\cos \theta}{2\pi}. \tag{28}\]

See Appendix B.1 for a derivation of these quantities. Note that Cho and Saul (2009) have previously discovered the formulae for \(J_{0,0}\) and \(J_{1,1}\) by use of a completely different contour-integral based method.

The combination of Proposition 1 and 2 make it possible to practically calculate any value of \(J_{a,b}\) when \(a,b\) are not too large. However, by serendipity, we are able to find remarkable explicit formulas for \(J_{a,b}\), which we report below.
**Proposition 3** (Explicit Formulas for \(J_{a,0}(\theta)\), \(J_{a,1}(\theta)\)): _Let \(a\geq 2\). Then, \(J_{a,0}\) and \(J_{a,1}\) are explicitly given by the following:_

\[J_{a,0}(\theta)=(a-1)!\,!\left(J_{a\bmod 2,0}+\frac{\cos\theta}{c_{a\bmod 2}} \sum_{\begin{subarray}{c}i\neq a(\bmod 2)\\ 0<i<a\end{subarray}}\frac{(i-1)!\,!}{i!\,!}\sin^{i}\theta\right),\]

_where \(c_{0}=2\pi\), \(c_{1}=2\sqrt{2\pi}\). We can then use the explicit formula for \(J_{a,0}\) in the formula for \(J_{a,1}\):_

\[J_{a,1}(\theta)=(a-1)!\,!\left(J_{a\bmod 2,1}+\cos\theta\sum_{ \begin{subarray}{c}i\neq a(\bmod 2)\\ 0<i<a\end{subarray}}\frac{J_{i,0}(\theta)}{i!\,!}\right),\]

_where an explicit formula for the first term (either \(J_{1,0}\) or \(J_{1,1}\) depending on the parity of \(a\)) is given in Proposition 2._

We can finally express \(J_{a,b}\) as a linear combination of \(J_{0,n}\) and \(J_{1,n}\), as follows. (In light of the previous explicit formulas, this is an explicit formula for \(J_{a,b}\).) It turns out that the coefficients are given in terms of two special numbers \(P(a,b)\) and \(Q(a,b)\) which are known as the Bessel numbers.

**Definition 3** (Bessel numbers): _The numbers \(P(a,b)\) and \(Q(a,b)\) are defined as follows,_

\[P(a,b) =\begin{cases}\frac{a!}{b!\left(\frac{a-b}{2}\right)^{2}\frac{a- b}{2}},&a\geq b,\;a\equiv b\;(\bmod 2)\\ 0,&\text{otherwise}\end{cases} \tag{29}\] \[Q(a,b) =\begin{cases}\frac{(\frac{a+b}{2})!}{b!}2\frac{b-a}{2}\sum_{i=0 }^{\frac{a-b}{2}}\binom{a+1}{i},&a\geq b,\;a\equiv b\;(\bmod 2)\\ 0,&\text{otherwise}\end{cases}. \tag{30}\]

\(P(a,b)\) represents a family of numbers known as the Bessel numbers of the second kind (Cheon et al. (2013)), and \(Q(a,b)\) comes from a closely related family of numbers (Kreinin (2016)). Using these, we can express \(J_{a,b}\) as follows.

**Theorem 4** (Explicit Formula for \(J_{a,b}(\theta)\)): _Let \(b\geq 2,a\geq 1,\;b\geq a\). Then, we have the following formula for \(J_{a,b}(\theta)\) in terms of \(J_{0,n}\) and \(J_{1,n}\)_

\[J_{a,b} =\sum_{\begin{subarray}{c}i\equiv 0(\bmod 2)\\ 0<i\leq a\end{subarray}}(b)_{a-i}(\cos\theta)^{a-i}\left(P(a,a-i)-Q(a-1,a-1-i )\right)J_{0,b-a+i}\] \[+\sum_{\begin{subarray}{c}i\equiv 1(\bmod 2)\\ 0<i\leq a\end{subarray}}(b)_{a-i}(\cos\theta)^{a-i}Q(a-1,a-i)J_{1,b-a+i}.\]
**Remark 5**: _Since \(J_{1,n}\) is also given in terms of \(J_{0,n}\), one may further simplify the formula for \(J_{a,b}\) to be in terms of only \(J_{0,n}\), \(J_{0,0}\) and \(J_{0,1}\). This substitution yields the following formula. For notational convenience, we will let \(\delta:=b-a\),_

\[J_{a,b} =\sum_{\begin{subarray}{c}i\equiv 0(\bmod 2)\\ 0<i\leq a\end{subarray}}(b)_{a-i}(\cos\theta)^{a-i}(P(a,a-i)-Q(a-1,a-1-i))J_{0,\delta+i}\] \[+\sum_{\begin{subarray}{c}i\equiv 1(\bmod 2)\\ 0<i\leq a\end{subarray}}(b)_{a-i}(\cos\theta)^{a-i}Q(a-1,a-i)(\delta+i-1)! \,J_{(\delta+1)\bmod 2,1}\] \[+\cos\theta\sum_{\begin{subarray}{c}i\equiv 1(\bmod 2)\\ 0<i\leq a\end{subarray}}\sum_{\begin{subarray}{c}j\equiv\delta(\bmod 2)\\ 0<j<\delta+i\end{subarray}}(b)_{a-i}(\cos\theta)^{a-i}Q(a-1,a-i)\frac{(\delta+i- 1)!}{j!\,!}J_{0,j}.\]

### Gaussian Integration-by-Parts Formulae

In this section, we state two important formulas that together give us the tools for computing the expectations that appear in \(J_{a,b}\).

**Fact 1** (Gaussian Integration by Parts): _Let \(G\sim\mathcal{N}(0,1)\) be a Gaussian variable and \(f:\mathbb{R}\to\mathbb{R}\) be a differentiable function. Then,_

\[\mathbf{E}[Gf(G)]=\mathbf{E}[f^{\prime}(G)]. \tag{31}\]

Using this type of Gaussian integration by parts formula, we can generalize the expected value of Gaussians to derivatives of functions which are not necessarily differentiable. For example the indicator function \(1\{x>a\}\) is not differentiable, but for the purposes of computing Gaussian expectation, we can use the following integration formula.

**Fact 2** (Gaussian expectations involving \(1^{\prime}\{x>a\}\)): _Let \(G\) be a Gaussian variable and \(a\in\mathbb{R}\). Let \(f:\mathbb{R}\to\mathbb{R}\) such that \(\lim_{g\to\infty}f(g)e^{\frac{-g^{2}}{2}}=0\). Then, using the Gaussian integration by parts formula to assign a meaning to expectations involving the "derivative of the indicator function", \(1^{\prime}\{x>a\}\), we have_

\[\mathbf{E}[1^{\prime}\{G>a\}f(G)]=f(a)\frac{e^{\frac{-g^{2}}{2}}}{\sqrt{2\pi}}. \tag{32}\]

**Remark 6**: _The purpose of assigning a value to the expectation (32) is to allow one to compute "honest" expectations of the form (31) when \(f(x)\) involves \(1\{x>0\}\); see Lemma 7 for an illustrative example. The final result does not require interpreting "\(1^{\prime}\{x>a\}\)"; this is only a useful intermediate step in the sequence of calculations leading to the final result._

_The formula can also be understood or proven in a number of different alternative ways. One is simply to say that \(1^{\prime}\{x>a\}=\delta\{x=a\}\) is a "Dirac delta function" at \(x=a\). A more rigorous way would be to take any differentiable family of functions \(1_{\epsilon}\{x>a\}\) which suitably converge to \(1\{x>a\}\) as \(\epsilon\to 0\) and then interpret the result as the limit of the expectation \(\lim_{\epsilon\to 0}\mathbf{E}[1^{\prime}_{\epsilon}\{G>a\}f(G)]\)._
**Proof** [Of Fact 2] Applying integration by parts, we formally have

\[\mathbf{E}[1^{\prime}\{G>a\}f(G)] =\int\limits_{-\infty}^{\infty}1^{\prime}\{g>a\}f(g)\frac{e^{\frac {-g^{2}}{2}}}{\sqrt{2\pi}}dg\] \[=\left[1\{g>a\}f(g)\frac{e^{\frac{-g^{2}}{2}}}{\sqrt{2\pi}} \right]_{-\infty}^{\infty}-\int\limits_{-\infty}^{\infty}1\{g>a\}\frac{d}{dg} \left(f(g)\frac{e^{\frac{-g^{2}}{2}}}{\sqrt{2\pi}}\right)dg.\]

Note that the first term is \(0\) by the hypothesis \(\lim\limits_{g\to\infty}f(g)e^{\frac{-g^{2}}{2}}=0\), and we have then

\[\mathbf{E}[1^{\prime}\{G>a\}f(G)]=-\int\limits_{a}^{\infty}\frac{d}{dg}\left( f(g)\frac{e^{\frac{-g^{2}}{2}}}{\sqrt{2\pi}}\right)dg=-\left[f(g)\frac{e^{ \frac{-g^{2}}{2}}}{\sqrt{2\pi}}\right]_{a}^{\infty}=0+f(a)\frac{e^{\frac{-g^{ 2}}{2}}}{\sqrt{2\pi}},\]

where we have used the hypothesis on \(f\) once again.

The two facts about Gaussian integration by parts can be combined to create recurrence relations for expectations involving \(1\{G>a\}\). A simple example is the following lemma, which we will also use later in our derivation. The proof strategy of this lemma is a microcosm of the proof strategy we use to compute \(J_{a,b}\) in general, namely to use Gaussian integration by parts to derive a recurrence relation and initial condition, and then solve.

**Lemma 7** (Moments of \(\varphi(G)\)): _For \(k\geq 0\), we have_

\[\mathbf{E}[\varphi(G)^{k}]=\mathbf{E}[G^{k}1\{G>0\}]=\begin{cases}\frac{(k-1)!! }{2}&\text{$k$ is even}\\ \frac{(k-1)!!}{\sqrt{2\pi}}&\text{$k$ is odd}\end{cases}=\sqrt{2\pi}\frac{(k-1)!! }{c_{k-1\bmod 2}},\]

_where \(c_{0}=2\pi\) and \(c_{1}=2\sqrt{2\pi}\)._

**Proof** We prove this for even and odd \(k\) separately by induction on \(k\). The base case for \(k=0\) is trivial since \((0-1)!!=1\) is the empty product. The base case \(k=1\) follows by first applying (31) with \(f(x)=1\{x>0\}\) and then applying (32) with \(f(x)\equiv 1\),

\[\mathbf{E}[\varphi(G)]=\mathbf{E}[G1\{G>0\}]=\mathbf{E}[1^{\prime}\{G>0\}]= \sqrt{2\pi}^{-1}.\]

Now, to see the induction, we apply (31) with \(f(x)=x^{k-1}1\{x>0\}\), \(k\geq 2\). Due to the product rule, there are two terms in the derivative,

\[\mathbf{E}[\varphi(G)^{k}] =\mathbf{E}[G\cdot G^{k-1}1\{G>0\}] \tag{33}\] \[=(k-1)\mathbf{E}[G^{k-2}1\{G>0\}]+\mathbf{E}[G^{k-1}1^{\prime}\{G >0\}]\] (34) \[=(k-1)\mathbf{E}[\varphi(G)^{k-2}]+0, \tag{35}\]

where we have recognized that the second term is \(0\) by application of (32) with \(f(x)=x^{k-1}\) which has \(f(0)=0\). The recurrence \(\mathbf{E}[\varphi(G)^{k}]=(k-2)\mathbf{E}[\varphi(G)^{k-2}]\) along with initial condition leads to the stated result by induction.
### Recursive Formulae for \(J_{a,b}(\theta)\) - Proof of Proposition 1

**Proof** [Of Proposition 1] To find a recursive formula for \(J_{a,0},a\geq 2\), we apply the Gaussian integration by parts formula (31) to \(f(x)=x^{a-1}1\{x>0\}1\{\cos\theta x+W\sin\theta>0\}\) to evaluate the expected value over \(G\) first. When applying product rule there are three terms

\[J_{a,0}= \mathbf{E}[G\cdot G^{a-1}1\{G>0\}1\{G\cos\theta+W\sin\theta>0\}] \tag{36}\] \[= (a-1)\mathbf{E}[G^{a-2}1\{G>0\}1\{G\cos\theta+W\sin\theta>0\}]\] \[+\mathbf{E}[G^{a-1}1\{G>0\}1^{\prime}\{G\cos\theta+W\sin\theta>0 \}]\cos\theta\] \[+\mathbf{E}[G^{a-1}1^{\prime}\{G>0\}1\{G\cos\theta+W\sin\theta>0 \}].\]

The first term is simply \((a-1)J_{a-2,0}\). The last two terms can now be evaluated with the help of (32). The last term of (36) is (32) with the function \(f(x)=x^{a-1}1\{x\cos\theta+W\sin\theta>0\}\) which has \(f(0)=0\) for \(a\geq 2\). Therefore, this term simply vanishes.

To evaluate the middle term of (36), we introduce a change of variables to express \(G\cos\theta+W\sin\theta\) in terms of two other independent Gaussian variables \(Z,W\sim\mathcal{N}(0,1)\)

\[Z =G\cos\theta+W\sin\theta G=Z\cos\theta+Y\sin\theta \tag{37}\] \[Y =G\sin\theta-W\cos\theta W=Z\sin\theta-Y\cos\theta,\]

where \(Y,Z\) iid \(\mathcal{N}(0,1)\). Under this change of variables, \(J_{a,0},\ a\geq 2\) is setup to apply (32) with \(f(x)=1\{x\cos\theta+Y\sin\theta\}^{a-1}1\{x\cos\theta+Y\sin\theta>0\}\)

\[J_{a,0} =(a-1)J_{a-2,0}+\mathbf{E}[G^{a-1}1\{G>0\}1^{\prime}\{G\cos \theta+W\sin\theta>0\}]\cos\theta\] \[=(a-1)J_{a-2,0}+\mathbf{E}[(Z\cos\theta+Y\sin\theta)^{a-1}1\{Z \cos\theta+Y\sin\theta>0\}1^{\prime}\{Z>0\}]\cos\theta\] \[=(a-1)J_{a-2,0}+\mathbf{E}[(0+Y\sin\theta)^{a-1}1\{0+Y\sin\theta> 0\}]\frac{1}{\sqrt{2\pi}}\cos\theta\] \[=(a-1)J_{a-2,0}+\frac{\sin^{a-1}\theta\cos\theta}{c_{a\bmod 2}}( a-2)!!\,,\]

where we have applied Lemma 7 to evaluate the last expectation.

A similar argument is used to find the recursive formula for \(J_{a,b},\ a\geq 2,b\geq 1\), by using (31) with the function \(f(x)=x^{a-1}(x\cos\theta+W\sin\theta)^{b}1\{x>0\}1\{x\cos\theta+W\sin\theta>0\}\). There are 4 terms in the product rule derivative. Fortunately in this case, the last two terms are simply zero by application of (32) since the expressions vanish when \(G=0\), so we get

\[J_{a,b} =\mathbf{E}[G\cdot G^{a-1}(G\cos\theta+W\sin\theta)^{b}1\{G>0\}1 \{G\cos\theta+W\sin\theta>0\}]\] \[=\mathbf{E}[(a-1)G^{a-2}(G\cos\theta+W\sin\theta)^{b}1\{G>0\}1\{G \cos\theta+W\sin\theta>0\}]\] \[\qquad+\mathbf{E}[G^{a-1}b\cos\theta(G\cos\theta+W\sin\theta)^{b- 1}1\{G>0\}1\{G\cos\theta+W\sin\theta>0\}]\] \[\qquad+\mathbf{E}[G^{a-1}(G\cos\theta+W\sin\theta)^{b}1^{\prime} \{G>0\}1\{G\cos\theta+W\sin\theta>0\}]\] \[\qquad+\mathbf{E}[G^{a-1}(G\cos\theta+W\sin\theta)^{b}1\{G>0\}1^{ \prime}\{G\cos\theta+W\sin\theta>0\}\cos\theta]\] \[=(a-1)J_{a-2,b}+b\cos\theta J_{a-1,b-1}+0+0,\]

as desired.
### Solving the Recurrence to get an Explicit Formula for \(J_{a,b}(\theta)\) - Proof of Theorem 4

Solving the recurrence for the sequences \(J_{a,0}\) and \(J_{a,1}\) to get the claimed explicit formula for \(J_{a,0}\) is a simple induction proof. We defer these to Appendix B.2. More difficult and interesting is the 2D array \(J_{a,b}\). To solve the recurrence

\[J_{a,b}=(a-1)J_{a-2,b}+b\cos\theta J_{a-1,b-1},\quad a\geq 2,b\geq 1, \tag{38}\]

we will apply the recursion repeatedly until \(J_{a,b}\) can be expressed as a linear combination of \(J_{0,n}\) and \(J_{1,n}\) terms for which we already have an explicit formula developed. To determine the coefficients in front of \(J_{0,n}\) and \(J_{1,n}\), we take a combinatorial approach by thinking of the recurrence relation as a weighted directed graph as defined below.

**Definition 8** (Viewing a recursion as a directed weighted graph): _We can view the recurrence relation for \(J_{a,b}\) as a weighted directed graph on the vertex set \((a,b)\in\mathbb{Z}^{2}\) where vertices represent the values of \(J_{a,b}\) and directed edges capture how values of \(J_{a,b}\) are connected through the recurrence relation. To be precise, the graph edges and edge weights \(w_{e}\) are defined so that the recursion (38) for \(J_{a,b}\) can be expressed in the graph as a sum over incoming edges,_

\[J_{a,b}=\sum_{e:(a^{\prime},b^{\prime})\to(a,b)}w_{e}^{J}J_{a^{\prime},b^{ \prime}}, \tag{39}\]

_where the sum is over the edges \(e\) with weight \(w_{e}^{J}\) incoming to the vertex \((a,b)\). An example of the graph to calculate \(J_{6,8}\) is illustrated in Figure 4._

_By repeatedly applying the recursion, \(J_{a,b}\) can be expressed as a linear combination of the values at the source vertices of the graph (i.e. those with no incoming edges). For the recurrence \(J_{a,b}\), the source vertices are \(J_{0,n}\) and \(J_{1,n}\). The coefficients in front of each source is simply the weighed sum over all paths from the source to the node, namely_

\[J_{a,b}=\sum_{\text{source vertices }v}W_{v\to(a,b)}^{J}J_{v} =\sum_{n\geq 0}W_{(0,n)\to(a,b)}^{J}J_{0,n}+\sum_{n\geq 0}W_{(1,n)\to(a,b)}^{J}J_{1,n}, \tag{40}\] \[W_{(a^{\prime},b^{\prime})\to(a,b)}^{J} :=\sum_{\pi:(a^{\prime},b^{\prime})\to(a,b)}\prod_{e\in\pi}w_{e}^ {J}, \tag{41}\]

_where the sum is over all paths \(\pi\) from the vertex \((a^{\prime},b^{\prime})\) to the vertex \((a,b)\) in the \(J\) graph._

In light of (40), to prove Theorem 4, we have only to calculate the weighted sum of path \(W_{(0,n)\to(a,b)}^{J}\) and \(W_{(1,n)\to(a,b)}^{J}\). These weighted sums turn out to be given in terms of the \(P\) and \(Q\) numbers which were defined in Definition 3.

**Proposition 4** (Weighted sums of paths for \(J\)): _In the graph for \(J\), we have the following formulas for the sum over weighted paths \(W^{J}\) defined in (41),_

\[W_{(0,n)\to(a,b)}^{J} =(b)_{b-n}(\cos\theta)^{b-n}(P(a,b-n)-Q(a-1,b-n-1)), \tag{42}\] \[W_{(1,n)\to(a,b)}^{J} =(b)_{b-n}(\cos\theta)^{b-n}Q(a-1,b-n). \tag{43}\]
To prove this Proposition 4, we first create a simpler recursion, \(J^{*}\), which we solve first and then slightly modify the solution to get the solution for \(J\).

**Lemma 9** (Weighted sums of paths for \(J^{*}\)): _Let \(J^{*}_{a,b}\) be defined to be the recursion:_

\[J^{*}_{a,b} :=(a-1)J^{*}_{a-2,b}+1J^{*}_{a-1,b-1}\text{ for }2\leq a\leq b \tag{44}\] \[J^{*}_{1,b} :=0+1J^{*}_{0,b-1}\text{ for }1\leq b. \tag{45}\]

_Thinking of this recursion as a graph as in Definition 8 (see Figure 4 for an illustration), we have that the sum of weighted paths \(W^{J^{*}}\) defined analogously to those in (41), are given by \(P\) and \(Q\) numbers, namely_

\[W^{J^{*}}_{(0,n)\to(a,b)} =P(a,b-n) \tag{46}\] \[W^{J^{*}}_{(1,n)\to(a,b)} =Q(a-1,b-n). \tag{47}\]

The connection between the \(J^{*}\) and the \(P\), \(Q\) numbers is through the following recursion for the \(P\), \(Q\) numbers.

**Lemma 10**: _(Recursion for \(P\) and \(Q\) numbers) The \(P\) numbers, defined in Definition 3, satisfy \(P(0,0)=1\), \(P(n,n)=P(n-1,n-1)\) for \(n\geq 1\), and the recursion_

\[P(a,b)=(\texttt{a}-1)\cdot P(a-2,b)+1\cdot P(a-1,b-1),\ \text{ for }\ a\geq 2,\ 0\leq b\leq a-2,\]

_under the convention that \(P(a,-1)=0\). The \(Q\) numbers satisfy the same recursion as the \(P\) numbers, with a coefficient of a rather than \((\texttt{a}-1)\)._

The proof of Lemma 10 is an easy consequence of known results from Kreinin (2016) and is deferred to Appendix C.

Figure 4: The graph associated with the recursions for \(J\) in (38) (left) and \(J^{*}\) in (45) (right). The graph is defined so that the recursion is given by a sum of incoming edges as in (39). The edges are color coded red and blue to match the coefficients in the recursion.
**Proof** [Of Lemma 9] Using the same idea of recursions expressed as graphs as in Definition 8, the recursion from Lemma 10 means that \(P\) and \(Q\) can be expressed as weighted directed graphs. These are displayed in Figure 5. Since the \(P\) and \(Q\) graphs have only one single unit valued source vertex at \((0,0)\), (40) shows that the \(P\) and \(Q\) numbers are actually themselves equal to sums over weighted paths in their respective graphs

\[P(a,b)=W^{P}_{(0,0)\rightarrow(a,b)},\quad Q(a,b)=W^{Q}_{(0,0)\rightarrow(a,b)}.\]

Therefore the statement of the lemma is that sum over weighted paths in the \(J^{*}\) graph are the same as other sums over weighted paths in the \(P\) graph/\(Q\) graphs,

\[W^{J^{*}}_{(0,n)\rightarrow(a,b)}=W^{P}_{(0,0)\rightarrow(a,b-n)},\quad W^{ J^{*}}_{(1,n)\rightarrow(a,b)}=W^{Q}_{(0,0)\rightarrow(a-1,b-n)}. \tag{48}\]

The fact that these are equal is demonstrated by establishing a simple bijection between weighted paths in the \(P\) graph/\(Q\) graph, and weighted paths in the \(J^{*}\) graph.

For example, in Figure 5, there is a bijection between the weighted paths in the \(P\) graph which connect \(P(0,0)\) to \(P(6,2)\), to the paths which connect \(J^{*}_{6,8}\) to \(J^{*}_{0,6}\) in the \(J^{*}\) graph. The bijection is simply to _flip_ any path in the \(P\)-graph by rotating it by \(180^{\circ}\) to get a valid path in the \(J^{*}\)-graph. Moreover, the edge weights for \(J^{*}\) and \(P\) are precisely set up so that under this bijection, the paths will have the same set of weighted edges in the same order. A full, more detailed, explanation of this bijection is given in Appendix B.3. This argument shows that \(W^{J^{*}}_{(0,n)\rightarrow(a,b)}=P(a,b-n)\) as desired.

Figure 5: Graphs associated with the recursions for the \(P\) numbers (left), \(Q\) numbers (middle), and \(J^{*}\) (right). The weighted edges indicate the coefficients in the recursions for \(P,Q,J^{*}\) respectively. The diagrams are lined up so that the sum of weighted paths in from \(J^{*}_{6,8}\) can be directly read from the \(P\) and \(Q\) entries in the same location. By reading from the bottom displayed row of \(P\) we see that the weighted sum over paths \(W^{J^{*}}_{(0,n)\rightarrow(6,8)}\) are \(15\), \(45\), \(15\) and \(1\) for \(n=8,6,4\) and \(2\) respectively. (Since these are the source vertices, this shows that \(J^{*}_{6,8}=15J^{*}_{0,8}+45J^{*}_{0,6}+15J^{*}_{0,4}+1J^{*}_{0,2}\).) From the bottom displayed row of \(Q\) we see that the values for sums of weighted paths from vertices \(W^{J^{*}}_{(1,n)\rightarrow(6,8)}\) are \(33\), \(14\) and \(1\) for \(n=7,5\) and \(3\) respectively.
The \(P\) numbers do not apply for paths between \(J^{*}_{1,n}\) and \(J^{*}_{a,b}\) because we are starting one row higher so the first vertical upward edge is weight 2. In this case, there is a bijection to the \(Q\)-graph after flipping the path. For any path which runs from a node in row 1 to the top left corner of the \(J^{*}\)-graph, we can find the same "flipped" path in the graph of the \(Q\), running from the top left entry to the corresponding node in row \(a-1\). (The bijection is explained in detail in Appendix B.3.) Hence \(W^{J^{*}}_{(1,n)\to(a,b)}=Q(a-1,b-n)\) as desired.

Having solved for \(J^{*}\) in terms of \(P\) and \(Q\), it remains to translate these into the weights for \(J\) to obtain Proposition 4.

**Proof** [Of Proposition 4]

The proof follows by relating the weighted sum of paths for \(J\) in terms of \(J^{*}\) and then applying the result of Lemma 9. There are two differences between the formula for \(J_{a,b}\) compared to \(J^{*}_{a,b}\), which can both be seen in Figure 4. We handle both differences as follows:

_Difference #1:_\(J\) has a weight of \(b\cos\theta\) on the blue diagonal edges \((a,b)\to(a+1,b+1)\) vs \(J^{*}\) has a weight of 1.

This difference is handled by the following observation: any path from \((a,b)\to(a^{\prime},b^{\prime})\) in the graph goes through each column between \(b\) and \(b^{\prime}\) exactly once. This means that the contribution of the edge weights from these edges do not depend on the details of which path was taken, only the starting and ending points. They always contribute the same factor, \((b)_{b^{\prime}-b}(\cos\theta)^{b^{\prime}-b}\). (here \((b)_{k}=b(b-1)\cdots(b-k+1)\) is the falling factorial with \(k\) terms). This argument shows that the weighted sum of paths in \(J\) and \(J^{*}\) are related by

\[W^{J}_{(a,b)\to(a^{\prime},b^{\prime})}=(b)_{b^{\prime}-b}(\cos\theta)^{b^{ \prime}-b}W^{J^{*}}_{(a,b)\to(a^{\prime},b^{\prime})}. \tag{49}\]

By the result of Lemma 9, this shows that \(W^{J}_{(1,n)\to(a,b)}=(b)_{b-n}(\cos\theta)^{b-n}Q(a-1,b-n)\) as desired.

Equation (49) holds for all paths with starting point \(a\geq 1\). When \(a=0\), there is one additional difference between \(J\) and \(J^{*}\) which is accounted for below.

_Difference #2:_\(J_{0,n}\) has _no_ diagonal edge vs \(J^{*}_{0,n}\) has a diagonal blue edge of weight 1.

Because of this "missing edge", the only choice in \(J\) for paths starting from \((0,n)\) is to first go vertically up by 2 units to \((2,n)\). Hence \(W^{J}_{(0,n)\to(a,b)}=W^{J}_{(2,n)\to(a,b)}\). To evaluate this, we use the decomposition of paths in \(J^{*}\) by what their first step is, either a diagonal blue step or a red vertical up step, to see that

\[W^{J^{*}}_{(0,n)\to(a,b)} =W^{J^{*}}_{(1,n+1)\to(a,b)}+W^{J^{*}}_{(2,n)\to(a,b)}, \tag{50}\] \[\implies W^{J^{*}}_{(2,n)\to(a,b)} =W^{J^{*}}_{(0,n)\to(a,b)}-W^{J^{*}}_{(1,n+1)\to(a,b)}\] (51) \[=P(a,b-n)-Q(a-1,b-n-1), \tag{52}\]

by the result of Lemma 9. By applying now (49) to relate \(J\) and \(J^{*}\), we obtain \(W^{J}_{(0,n)\to(a,b)}=W^{J}_{(2,n)\to(a,b)}=(b)_{b-n}(\cos\theta)^{b-n}(P(a,b -n)-Q(a-1,b-n-1))\) as desired.

**Proof** [Of Theorem 4] The formula is immediate from (40), which writes \(J_{a,b}\) as a linear combination of \(J_{0,n}\) and \(J_{1,n}\), and Proposition 4 which gives the the coefficients.
## Appendix A.

### Expected Value Approximation

**Lemma 11**: _Both the random variables \(X=R^{\ell+1}\) and \(X=R^{\ell+1}\sin^{2}(\theta^{\ell})\) satisfy_

\[\mathbf{E}[\ln(X)]=\ln(\mathbf{E}[X])-\frac{\mathbf{Var}[X]}{2 \mathbf{E}[X]^{2}}+\mathcal{O}(n_{\ell}^{-2}). \tag{53}\]

**Proof** First note that by the properties of the logarithm, we have

\[\ln(X)=\ln\left(\mathbf{E}[X]\left(\frac{\mathbf{E}[X]+(X-\mathbf{E}[X])}{ \mathbf{E}[X]}\right)\right)=\ln(\mathbf{E}[X])+\ln\left(1+\frac{X-\mathbf{E} [X]}{\mathbf{E}[X]}\right). \tag{54}\]

We can now apply the Taylor series \(\ln(1+x)=x-\frac{x^{2}}{2}+\epsilon_{2}(x)\), where \(\epsilon_{2}(x)\) is the Taylor series remainder and satisfies \(\epsilon_{2}(x)=\mathcal{O}(x^{3})\). Hence

\[\ln(X)=\ln(\mathbf{E}[X])+\frac{X-\mathbf{E}[X]}{\mathbf{E}[X]}- \frac{(X-\mathbf{E}[X])^{2}}{2\mathbf{E}[X]^{2}}+\epsilon_{2}\left(\frac{X- \mathbf{E}[X]}{\mathbf{E}[X]}\right).\]

Note that \(\mathbf{E}[X-\mathbf{E}[X]]=0\), and \(\mathbf{E}[(X-\mathbf{E}[X])^{2}]=\mathbf{Var}[X]\). Thus, if we take the expected value of our above approximation, we get the following:

\[\mathbf{E}[\ln(X)]=\ln(\mathbf{E}[X])-\frac{\mathbf{Var}[X]}{2 \mathbf{E}[X]^{2}}+\mathbf{E}\left[\epsilon_{2}\left(\frac{X-\mathbf{E}[X]}{ \mathbf{E}[X]}\right)\right].\]

By using bounds on the Taylor series error term \(\epsilon_{2}(x)=\mathcal{O}(x^{3})\), one can obtain bounds for this last error term. By (16, 17), both \(X=R^{\ell+1}\) and \(X=R^{\ell+1}\sin^{2}(\theta^{\ell+1})\) can be expressed as averages of the form

\[X=\frac{1}{n_{\ell}^{2}}\sum_{i,j}^{n_{\ell}}f(G_{i},\hat{G}_{j}). \tag{55}\]

From the bound on the 3rd moment in Lemma 14 it follows that \(\mathbf{E}[\epsilon_{2}(X-\mathbf{E}[X])]=\mathcal{O}(n_{\ell}^{-2})\), thus giving the desired result.

### Variance Approximation

**Lemma 12**: _Both the random variables \(X=R^{\ell+1}\) and \(X=R^{\ell+1}\sin^{2}(\theta^{\ell})\) satisfy_

\[\mathbf{Var}[\ln(X)]=\frac{\mathbf{Var}[X]}{\mathbf{E}[X]^{2}}+ \mathcal{O}(n_{\ell}^{-2}).\]

**Proof** Starting with (54), and using the first term of the Taylor series approximation for \(\ln(1+x)=x+\epsilon_{1}(x)\) now, we have that

\[\ln(X)=\ln(\mathbf{E}[X])+\frac{X-\mathbf{E}[X]}{\mathbf{E}[X]}+ \epsilon_{1}\left(\frac{X-\mathbf{E}[X]}{\mathbf{E}[X]}\right). \tag{56}\]
where \(\epsilon_{1}(x)\) is the Taylor error term and satisfies \(\epsilon_{1}(x)=\mathcal{O}(x^{2})\). Taking the variance of this, we arrive at an approximation of \(\mathbf{Var}[\ln(X)]\).

\[\mathbf{Var}[\ln(X)]=\mathbf{Var}\left[\ln(\mathbf{E}[X])+\frac{X-\mathbf{E}[X ]}{\mathbf{E}[X]}+\epsilon_{1}\left(\frac{X-\mathbf{E}[X]}{\mathbf{E}[X]} \right)\right]\]

\[=\mathbf{Var}\left[\frac{X-\mathbf{E}[X]}{\mathbf{E}[X]}\right]+\mathbf{Var} \left[\epsilon_{1}\left(\frac{X-\mathbf{E}[X]}{\mathbf{E}[X]}\right)\right]+ \mathbf{Cov}\left(\frac{X-\mathbf{E}[X]}{\mathbf{E}[X]},\epsilon_{1}\left( \frac{X-\mathbf{E}[X]}{\mathbf{E}[X]}\right)\right).\]

As with the expected value approximation, this approximation for variance is used twice, once for \(X=R^{\ell+1}\), and once for \(X=R^{\ell+1}\sin^{2}(\theta^{\ell+1})\) (see Section 2.2), both of which can can be expressed as a sum as in (55). Since \(\epsilon_{1}(x)=\mathcal{O}(x^{2})\), we have that the terms with \(\epsilon_{1}(x)\) are both \(\mathcal{O}(n_{\ell}^{-2})\) from Lemma 14. Simplifying the first term, \(\mathbf{Var}\left[\frac{X-\mathbf{E}[X]}{\mathbf{E}[X]}\right]=\frac{\mathbf{ Var}[X]}{\mathbf{E}[X]^{2}}\) gives the result of the Lemma. \(\blacksquare\)

### Covariance Approximation

**Lemma 13**: _Both the random variables \(X=R^{\ell+1}\) and \(X=R^{\ell+1}\sin^{2}(\theta^{\ell})\) satisfy_

\[\mathbf{Cov}(\ln(X),\ln(Y))=\frac{\mathbf{Cov}(X,Y)}{\mathbf{E}[X]\mathbf{E}[ Y]}+\mathcal{O}(n_{\ell}^{-2}).\]

**Proof** Using the approximation in (56) for \(\ln(X)\) and \(\ln(Y)\), we get the following expression for the covariance:

\[\mathbf{Cov}(\ln(X),\ln(Y))\] \[=\mathbf{Cov}\left(\frac{X}{\mathbf{E}[X]},\frac{Y}{\mathbf{E}[Y] }\right)+\mathbf{Cov}\left(\frac{X}{\mathbf{E}[X]},\epsilon_{1}\left(\frac{Y- \mathbf{E}[Y]}{\mathbf{E}[Y]}\right)\right)+\mathbf{Cov}\left(\epsilon_{1} \left(\frac{X-\mathbf{E}[X]}{\mathbf{E}[X]}\right),\frac{Y}{\mathbf{E}[Y]}\right)\] \[\quad+\mathbf{Cov}\left(\epsilon_{1}\left(\frac{X-\mathbf{E}[X]}{ \mathbf{E}[X]}\right),\epsilon_{1}\left(\frac{Y-\mathbf{E}[Y]}{\mathbf{E}[Y] }\right)\right).\]

We get the desired result from the fact that that the error term \(\epsilon_{1}(x)\) satisfies \(\epsilon_{1}(x)=\mathcal{O}(x^{2})\) and from our result in Lemma 16. \(\blacksquare\)

### Third and Fourth Moment Bound Lemma

**Lemma 14**: _Let \(G_{i},\hat{G}_{i}\), \(1\leq i\leq n\) be marginally \(\mathcal{N}(0,1)\) random variables with correlation \(\cos(\theta)\) and independent for different indices \(i\). Let \(A=\frac{1}{n^{2}}\sum_{i,j}^{n}f(G_{i},\hat{G}_{j})\) be the average over all \(n^{2}\) pairs of some function \(f:\mathbb{R}^{2}\rightarrow\mathbb{R}\) which has finite fourth moment, \(\mathbf{E}[f(G_{i},\hat{G}_{i})^{4}]<\infty\). Then the third and fourth central moment of \(A\) satisfy_

\[\mathbf{E}[(A-\mathbf{E}[A])^{3}] =\mathcal{O}(n^{-2}) \tag{57}\] \[\mathbf{E}[(A-\mathbf{E}[A])^{4}] =\mathcal{O}(n^{-2}). \tag{58}\]
**Proof** We begin by showing the third moment bound. First, we can express \({\bf E}[(A-{\bf E}[A])^{3}]\) as a sum in the following way:

\[A-{\bf E}[A]=\frac{1}{n^{2}}\sum_{i,j}^{n}\left(f(G_{i},\hat{G}_{j})-{\bf E}[f(G _{i},\hat{G}_{j})]\right)\]

\[\Longrightarrow\ {\bf E}\left[(A-{\bf E}[A])^{3}\right]=\frac{1}{n^{6}}\sum_{ i_{1},i_{2},i_{3}\atop j_{1},j_{2},j_{3}}^{n}{\bf E}\left[\prod_{k=1}^{3}\left(f(G_{i_ {k}},\hat{G}_{j_{k}})-{\bf E}[f(G_{i_{k}},\hat{G}_{j_{k}})]\right)\right]. \tag{59}\]

Note that many of these terms are mean zero. For example, for any configuration of the indices where there is no overlap between the indices \((i_{1},j_{1})\) and the other two index pairs, \(\{i_{1},j_{1}\}\cap\{i_{2},j_{2},i_{3},j_{3}\}=\emptyset\) we may use independence to observe that

\[{\bf E}\left[\prod_{k=1}^{3}\left(f(G_{i_{k}},\hat{G}_{j_{k}})-{ \bf E}[f(G_{i_{k}},\hat{G}_{j_{k}})]\right)\right]\] \[= {\bf E}\left[f(G_{i_{1}},\hat{G}_{j_{1}})-{\bf E}[f(G_{i_{1}}, \hat{G}_{j_{1}})]\right]{\bf E}\left[\prod_{k=2}^{3}\left(f(G_{i_{k}},\hat{G}_ {j_{k}})-{\bf E}[f(G_{i_{k}},\hat{G}_{j_{k}})]\right)\right]=0.\]

When this happens we say that \((i_{1},j_{1})\) is a "reducible point". Similarly, \((i_{2},j_{2})\) or \((i_{3},j_{3})\) can be reducible if they have no overlap with the other two index pairs. To control \({\bf E}\left[(A-{\bf E}[A])^{3}\right]\), it will suffice to enumerate the number of indices \(\{i_{1},j_{1},i_{2},j_{2},i_{3},j_{3}\}\) so that all three points \((i_{1},j_{1}),(i_{2},j_{2}),(i_{3},j_{3})\) are _not_ reducible. We call these "irreducible configurations".

We now observe that at least one of the points \((i_{1},j_{1})\),\((i_{2},j_{2})\) or \((i_{3},j_{3})\) is reducible whenever the number of unique numbers is \(\left|\bigcup_{k=1}^{3}\{i_{k},j_{k}\}\right|\geq 5\). This is because, by the pigeonhole principle, if there are no repeated or only one repeated number between 6 indices, then at least one of the 3 pairs \((i_{1},j_{1})\),\((i_{2},j_{2})\) or \((i_{3},j_{3})\) must consist of two unique numbers and therefore is a reducible point.

Since the only irreducible configurations can only have at most 4 unique numbers, the number of irreducible configurations is \({\cal O}(n^{4})\) as \(n\to\infty\). In fact, a detailed enumeration of the number of configurations reveals that the number of irreducible configurations is precisely

\[32(n)_{4}+68(n)_{3}+28(n)_{2}+1(n)_{1} \tag{60}\]

The leading term is 32 because there are 32 possible "patterns" for how the indices can be arranged to be both irreducible and contain exactly 4 unique numbers \(\left|\bigcup_{k=1}^{3}\{i_{k},j_{k}\}\right|=4\); these patterns are listed in Table 3. Each pattern contributes \((n)_{4}=n(n-1)(n-2)(n-3)\) possible index configurations by filling in the 4 unique numbers in all the possible ways. Similarly, there are respectively 68, 28, and 1 pattern(s) for irreducible configurations with 3,2 and 1 unique number(s) in them which each contribute \((n)_{3},(n)_{2}\) and \((n)_{1}\) configurations per pattern (Here \((n)_{k}\) denotes the falling factorial with \(k\) terms).

Since the number of irreducible configurations is \({\cal O}(n^{4})\), the normalization by \(n^{6}\) in (59) shows that \({\bf E}[(A-{\bf E}[A])^{3}]\) is \({\cal O}(n^{-2})\) as desired for the third moment.
The argument for the 4th moment is similar. We write \({\bf E}[(A-{\bf E}[A])^{4}]\) as a sum over \(i_{1},j_{1},i_{2},j_{2},i_{3},j_{3},i_{4},j_{4}\) and again enumerate irreducible configurations. In this case, once again by the pigeonhole principle any configuration with 7 or more unique points \(\left|\bigcup_{k=1}^{3}\{i_{k},j_{k}\}\right|\geq 7\) will be reducible. Since there are at most 6 unique numbers, there will be \({\cal O}(n^{6})\) irreducible configurations. A detailed enumeration of all the possible irreducible patterns and the number of unique elements in each yields that the number of irreducible configurations is precisely

\[48(n)_{6}+544(n)_{5}+1268(n)_{4}+844(n)_{3}+123(n)_{2}+1(n)_{1}\]

The normalization factor of \(n^{-8}\) then shows that \({\bf E}[(A-{\bf E}[A])^{4}]={\cal O}(n^{-2})\). \(\blacksquare\)

**Remark 15**: _A more detailed enumeration of the \(4\)-th moment actually shows that the dominant terms in the \(4\)th moment correspond to the terms in the \(2\)nd moment written twice, and asymptotically_

\[{\bf E}[(A-{\bf E}[A])^{4}]=3{\bf E}[(A-{\bf E}[A])^{2}]^{2}+{\cal O}(n^{-3}).\]

_Here \(3\) arises as the number of pair partitions of 4 items, and is related to the fact that \(3={\bf E}[G^{4}]\)_

**Lemma 16**: _Let \(G_{i},\hat{G}_{i},\,1\leq i\leq n\) be marginally \({\cal N}(0,1)\) random variables with correlation \(\cos(\theta)\) and independent for different indices \(i\). Let \(A_{1}=\frac{1}{n^{2}}\sum_{i,j}^{n}f_{1}(G_{i},\hat{G}_{j})\), and let \(A_{2}=\frac{1}{n^{2}}\sum_{i,j}^{n}f_{2}(G_{i},\hat{G}_{j})\), where \(f_{1},f_{2}:\mathbb{R}^{2}\rightarrow\mathbb{R}\) have finite fourth moments, \({\bf E}[f_{1}(G_{i},\hat{G}_{i})^{4}],\ {\bf E}[f_{2}(G_{i},\hat{G}_{i})^{4}]<\infty\). Then,_

\[{\bf E}[(A_{1}-{\bf E}[A_{1}])^{2}(A_{2}-{\bf E}[A_{2}])]={\cal O}\left(n^{-2} \right).\]

**Proof** We can express \({\bf E}[(A_{1}-{\bf E}[A_{2}])^{2}(A_{2}-{\bf E}[A_{2}])]\) using sums as follows:

\[{\bf E}[(A_{1}-{\bf E}[A_{1}])^{2}(A_{2}-{\bf E}[A_{2}])]\] \[=\frac{1}{n^{6}}\sum_{\begin{subarray}{c}i_{1},i_{2},i_{3}\\ j_{1},j_{2},j_{3}\end{subarray}}^{n}{\bf E}\left[\prod_{k=1}^{2}\left(f_{1}(G_{ i_{k}},\hat{G}_{j_{k}})-{\bf E}[f_{1}(G_{i_{k}},\hat{G}_{j_{k}})]\right)\left(f_{2}(G_ {i_{3}},\hat{G}_{j_{3}})-{\bf E}[f_{2}(G_{i_{3}},\hat{G}_{j_{3}})]\right) \right].\]

By the same argument as in Lemma 14, we can show that the number of nonzero terms in the above summation is \({\cal O}(n^{4})\) as \(n\rightarrow\infty\). Thus, we have that \({\bf E}[(A-{\bf E}[A])^{2}(A^{\prime}-{\bf E}[A^{\prime}])]={\cal O}(n^{-2})\). \(\blacksquare\)

### Derivation of Useful Identities - Equations (14, 15)

Let \(G\in\mathbb{R}^{n}\) be a Gaussian vector with iid entries \(G_{i}\sim{\cal N}(0,1)\). Then, by standard properties of Gaussians, the function \(f:\mathbb{R}^{n}\rightarrow\mathbb{R}\) given by \(f(x)=\langle G,x\rangle\) is a Gaussian random variable. Further, \(f(x)\) has the following properties:

* For all \(x\in\mathbb{R}^{n}\), \(f(x)\sim{\cal N}(0,\|x\|^{2})\)
* For any two vectors \(x_{\alpha},x_{\beta}\in\mathbb{R}^{n}\), the joint distribution of \(f(x_{\alpha}),f(x_{\beta})\) is jointly Gaussian with \[\begin{bmatrix}f(x_{\alpha})\\ f(x_{\beta})\end{bmatrix}\sim\mathcal{N}\left(0,\Sigma(x_{\alpha},x_{\beta}) \right),\qquad\Sigma(x_{\alpha},x_{\beta}):=\begin{bmatrix}\|x_{\alpha}\|^{2}& \langle x_{\alpha},x_{\beta}\rangle\\ \langle x_{\alpha},x_{\beta}\rangle&\|x_{\beta}\|^{2}\end{bmatrix}.\] \(\Sigma(x_{\alpha},x_{\beta})\) is sometimes called the \(2\times 2\) Gram matrix of the vectors \(x_{\alpha},x_{\beta}\).

In our setting of our fully connected neural network, any index \(i\in[n_{l+1}]\) in the vector of \(z^{\ell+1}\) is actually the inner product with the \(i\)-th row

\[z_{i}^{\ell+1}(x)=\sqrt{\frac{2}{n_{\ell}}}\langle W_{i,\cdot}^{\ell+1}, \varphi(z^{\ell}(x))\rangle.\]

Note that each row \(W_{i,\cdot}^{\ell+1}\) is a Gaussian vector, so the previous fact about Gaussians applies and we see that the entries of \(z^{\ell+1}\) are conditionally Gaussian given the value of the previous layer. By the previous Gaussian fact, we have that \(z_{i}^{\ell+1}(x_{\alpha})\), \(z_{i}^{\ell+1}(x_{\beta})\) are jointly Gaussian with

\[\begin{bmatrix}z_{i}^{\ell+1}(x_{\alpha})\\ z_{i}^{\ell+1}(x_{\beta})\end{bmatrix}\sim\mathcal{N}\left(0,\frac{2}{n_{\ell}} \begin{bmatrix}\|\varphi_{\alpha}^{\ell}\|^{2}&\langle\varphi_{\alpha}^{\ell},\varphi_{\beta}^{\ell}\rangle\\ \langle\varphi_{\alpha}^{\ell},\varphi_{\beta}^{\ell}\rangle&\|\varphi_{ \beta}^{\ell}\|^{2}\end{bmatrix}\right)=:\mathcal{N}\left(0,K^{\ell}\right),\]

where we use \(K^{\ell}\) to denote the \(2\times 2\) covariance matrix. \(K^{\ell}\) is precisely the \(2\times 2\) Gram matrix of the previous layer \(\varphi_{\alpha}^{\ell}\), \(\varphi_{\beta}^{\ell}\) scaled by \(2/n_{\ell}\) and its entries \(K_{\gamma\delta}^{\ell}\), for \(\gamma\in\{\alpha,\beta\},\delta\in\{\alpha,\beta\}\) are actually _averages_ of entries in the previous layer

\[K_{\gamma,\delta}^{\ell}:=\frac{2}{n_{\ell}}\langle\varphi_{\gamma}^{\ell}, \varphi_{\delta}^{\ell}\rangle=\frac{1}{n_{\ell}}\sum_{k=1}^{n_{\ell}}2\varphi (z_{k}^{\ell}(x_{i}))\varphi(z_{k}^{\ell}(x_{j})).\]

Moreover, in the weight matrix \(W^{\ell+1}\), the \(i^{th}\) and \(j^{th}\) rows (\(W_{i,\cdot}^{\ell+1}\) and \(W_{j,\cdot}^{\ell+1}\), respectively) are independent. Therefore, all entries of \(z^{\ell+1}\) are identically distributed and con

\begin{table}
\begin{tabular}{c c c c c} \((i_{1},j_{1})\) & \((i_{2},j_{2})\) & \((i_{3},j_{3})\) & \\ \hline \(\{(a,b),(b,a)\}\) & \(\times\) & \(\{(a,c),(c,a)\}\) & \(\times\) & \(\{(a,d),(d,a)\}\) & 8 patterns \\ \(\{(a,b),(b,a)\}\) & \(\times\) & \(\{(a,c),(c,a)\}\) & \(\times\) & \(\{(c,d),(d,c)\}\) & 8 patterns \\ \(\{(a,b),(b,a)\}\) & \(\times\) & \(\{(c,d),(d,c)\}\) & \(\times\) & \(\{(a,c),(c,a)\}\) & 8 patterns \\ \(\{(a,c),(c,a)\}\) & \(\times\) & \(\{(a,b),(b,a)\}\) & \(\times\) & \(\{(c,d),(c,b)\}\) & 8 patterns \\ \end{tabular}
\end{table}
Table 3: All 32 irreducible patterns using exactly 4 unique index values \(a,b,c,d\). For example the pattern \((i_{1},j_{1}),(i_{2},j_{2}),(i_{3},j_{3})=(a,b),(a,c),(a,d)\) represents all configurations where \(i_{1}=i_{2}=i_{3}\) and the \(j\)’s are all unique and different from \(i\). For each pattern, there are \((n)_{4}=n(n-1)(n-2)(n-3)\) configurations by filling in \(a,b,c,d\) with unique numbers in \([n]\). These are the dominant terms in (59).
ditionally independent given \(\varphi(z^{\ell})\). From this fact, we can equivalently write the entries explicitly as

\[z_{i}^{\ell+1}(x_{\alpha})=\sqrt{\frac{2}{n_{\ell}}}\|\varphi_{\alpha}^{\ell}\|G _{i},\quad z_{i}^{\ell+1}(x_{\beta})=\sqrt{\frac{2}{n_{\ell}}}\|\varphi_{\beta }^{\ell}\|\hat{G}_{i}, \tag{61}\]

where \(G_{i},\hat{G}_{i}\) are marginally \(\mathcal{N}(0,1)\) variables with covariance \(\mathbf{Cov}(G_{i},\hat{G}_{i})=\cos(\theta^{\ell})\) and independent for different indices. This formulation precisely ensures that the covariance structure for the entries is exactly what is specified by the covariance kernel \(K^{\ell}\).

With this representation of \(z_{i}^{\ell+1}(x_{\alpha})\) and \(z_{i}^{\ell+1}(x_{\beta})\), we can apply \(\varphi(\cdot)\) to each entry. By using the property of ReLU \(\varphi(\lambda x)=\lambda\varphi(x)\) for \(\lambda>0\) to factor out the norms, we obtain

\[\varphi(z_{i}^{\ell+1}(x_{\alpha}))=\sqrt{\frac{2}{n_{\ell}}}\|\varphi_{ \alpha}^{\ell}\|\varphi(G_{i}),\quad\varphi(z_{i}^{\ell+1}(x_{\beta}))=\sqrt{ \frac{2}{n_{\ell}}}\|\varphi_{\beta}^{\ell}\|\varphi(\hat{G}_{i}). \tag{62}\]

Taking the norm/inner product of the vector now yields (14-16) as desired.

### Cauchy-Binet and Determinant of the Gram Matrix - Equation (16)

To prove this identity, we begin with the fact that

\[\|\varphi_{\alpha}^{\ell+1}\|^{2}\|\varphi_{\beta}^{\ell+1}\|^{2}\sin^{2}( \theta^{\ell+1})=\det\begin{bmatrix}\|\varphi_{\alpha}^{\ell+1}\|^{2}&\langle \varphi_{\alpha}^{\ell+1},\varphi_{\beta}^{\ell+1}\rangle\\ \langle\varphi_{\alpha}^{\ell+1},\varphi_{\beta}^{\ell+1}\rangle&\|\varphi_{ \beta}^{\ell+1}\|^{2}\end{bmatrix}.\]

By the Cauchy-Binet identity, we can express the determinant as

\[\det\begin{bmatrix}\|\varphi_{\alpha}^{\ell+1}\|^{2}&\langle\varphi_{\alpha} ^{\ell+1},\varphi_{\beta}^{\ell+1}\rangle\\ \langle\varphi_{\alpha}^{\ell+1},\varphi_{\beta}^{\ell+1}\rangle&\|\varphi_{ \beta}^{\ell+1}\|^{2}\end{bmatrix}=\sum_{1\leq i<j\leq n_{\ell}}\left(\varphi_ {i;\alpha}^{\ell+1}\varphi_{j;\beta}^{\ell+1}-\varphi_{j;\alpha}^{\ell+1} \varphi_{i;\beta}^{\ell+1}\right)^{2}. \tag{63}\]

We can express layer \(\ell+1\) using the following conditioning on the previous layer

\[\varphi_{i;\alpha}^{\ell+1}=\sqrt{\frac{2}{n_{\ell}}}\|\varphi_{\alpha}^{\ell }\|\cdot\varphi(G_{i}),\quad\varphi_{i;\beta}^{\ell+1}=\sqrt{\frac{2}{n_{\ell}} }\|\varphi_{\beta}^{\ell}\|\cdot\varphi(\hat{G}_{i}).\]

Applying this to our expression in (63), and dividing both sides by \(\|\varphi_{\alpha}^{\ell}\|^{2}\|\varphi_{\beta}^{\ell}\|^{2}\), we get

\[\frac{\|\varphi_{\alpha}^{\ell+1}\|^{2}\|\varphi_{\beta}^{\ell+1}\|^{2}}{\| \varphi_{\alpha}^{\ell}\|^{2}\|\varphi_{\beta}^{\ell}\|^{2}}\sin^{2}(\theta^{ \ell+1})=\left(\frac{2}{n_{\ell}}\right)^{2}\sum_{1\leq i<j\leq n_{\ell}} \left(\varphi(G_{i})\varphi(\hat{G}_{j})-\varphi(G_{j})\varphi(\hat{G}_{i}) \right)^{2}.\]

Due to the fact that the summand is equal to \(0\) when \(i=j\), we can rewrite the sum as follows:

\[R^{\ell+1}\sin^{2}(\theta^{\ell+1})=\frac{1}{2}\left(\frac{2}{n_{\ell}}\right) ^{2}\sum_{i,j}^{n_{\ell}}\left(\varphi(G_{i})\varphi(\hat{G}_{j})-\varphi(G_{ j})\varphi(\hat{G}_{i})\right)^{2},\]

and we have the desired result.
### Expected Value calculations

In this section, we derive the formulas for \(\mathbf{E}\left[R^{\ell+1}\right]\), \(\mathbf{E}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right]\). We use \(J_{a,b}\) to represent \(J_{a,b}(\theta^{\ell})\). Note that \(\mathbf{E}[\varphi^{2}(G)]=\frac{1}{2}\), \(\mathbf{E}[\varphi^{4}(G)]=\frac{3}{2}\).

Calculation of \(\mathbf{E}\left[R^{\ell+1}\right]\) :

First, we apply the identity as in (17):

\[\mathbf{E}\left[R^{\ell+1}\right]=\left(\frac{2}{n_{\ell}}\right)^{2}\mathbf{ E}\left[\sum_{i,j=1}^{n_{\ell}}\varphi^{2}(G_{i})\varphi^{2}(\hat{G}_{j})\right].\]

Whenever \(i=j\), taking the expected value will give us a \(J_{2,2}\) term. When \(i\neq j\), the expected value of this term will be \(\mathbf{E}[\varphi^{2}(G)]^{2}=\frac{1}{4}\). Since \(i=j\) happens \(n_{\ell}\) times, and therefore \(i\neq j\) happens \(n_{\ell}^{2}-n_{\ell}\) times, we arrive at the following expression:

\[\mathbf{E}\left[R^{\ell+1}\right]=\left(\frac{2}{n_{\ell}}\right)^{2}\left(n _{\ell}J_{2,2}+(n_{\ell}^{2}-n_{\ell})\left(\frac{1}{4}\right)\right)=\frac{4J _{2,2}-1}{n_{\ell}}+1.\]

Calculation of \(\mathbf{E}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right]\) :

Applying the identity (16), we get

\[\mathbf{E}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right] =\frac{2}{n_{\ell}^{2}}\mathbf{E}\left[\sum_{i,j}^{n_{\ell}} \left(\varphi(G_{i})\varphi(\hat{G}_{j})-\varphi(G_{j})\varphi(\hat{G}_{i}) \right)^{2}\right]\] \[=\frac{2}{n_{\ell}^{2}}\mathbf{E}\left[\sum_{i,j}^{n_{\ell}} \left(\varphi^{2}(G_{i})\varphi^{2}(\hat{G}_{j})-2\varphi(G_{i})\varphi(\hat{ G}_{i})\varphi(G_{j})\varphi(\hat{G}_{j})+\varphi^{2}(G_{j})\varphi^{2}(\hat{G}_{i}) \right)\right].\]

In the case where \(i=j\), the expected value is equal to \(0\). Thus, we only need to consider the case where \(i\neq j\), which happens \(n_{\ell}^{2}-n_{\ell}\) times. When \(i\neq j\), the expectation of \(\varphi(G_{i})\varphi(\hat{G}_{i})\varphi(G_{j})\varphi(\hat{G}_{j})\) is \(J_{1,1}^{2}\), and the expectation of \(\varphi^{2}(G_{i})\varphi^{2}(\hat{G}_{j})\) is \(\frac{1}{4}\). All together, we have

\[\mathbf{E}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right]=\left(\frac{2}{n_{ \ell}^{2}}\right)(n_{\ell}^{2}-n_{\ell})\left(\frac{1}{4}-2J_{1,1}^{2}+\frac{ 1}{4}\right)=\frac{(n_{\ell}-1)(1-4J_{1,1}^{2})}{n_{\ell}}.\]

### Variance and Covariance Calculations

In this section, \(\mathbf{Var}\left[R^{\ell+1}\right]\), \(\mathbf{Var}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right]\), and \(\mathbf{Cov}\left(R^{\ell+1}\sin^{2}(\theta^{\ell+1}),R^{\ell+1}\right)\) are evaluated. We use \(J_{a,b}\) to represent \(J_{a,b}(\theta^{\ell})\). Note that \(\mathbf{E}[\varphi^{2}(G)]=\frac{1}{2}\), \(\mathbf{E}[\varphi^{4}(G)]=\frac{3}{2}\). We will see that there are simple functions \(f_{1},f_{2}:\mathbb{R}^{2}\rightarrow\mathbb{R}\) so that all of the variance and covariance calculations can be expressed as sums over \(i_{1},j_{1},i_{2},j_{2}\) of the form

\[\frac{1}{n_{\ell}^{4}}\sum_{\begin{subarray}{c}i_{1},j_{1}\\ i_{2},j_{2}\end{subarray}}\mathbf{Cov}\left[f_{1}(G_{i_{1}},\hat{G}_{j_{1}}), f_{2}(G_{i_{2}},\hat{G}_{j_{2}})\right] \tag{64}\] \[= \frac{1}{n_{\ell}^{4}}\sum_{\begin{subarray}{c}i_{1},j_{1}\\ i_{2},j_{2}\end{subarray}}\left(\mathbf{E}\left[f_{1}(G_{i_{1}},\hat{G}_{j_{1}}) f_{2}(G_{i_{2}},\hat{G}_{j_{2}})\right]-\mathbf{E}\left[f_{1}(G_{i_{1}},\hat{G}_{j_{1}}) \right]\mathbf{E}\left[f_{2}(G_{i_{2}},\hat{G}_{j_{2}})\right]\right), \tag{65}\]
where the sum goes over index configuration \((i_{1},j_{1}),(i_{2},j_{2})\in[n_{\ell}]^{4}\). We will use this form to organize our calculations of the variance and covariance formulas. The strategy is to evaluate each term in the sum (64) individually.

Since the random variables \(\{G_{i},\hat{G}_{i}\}_{i=1}^{n}\) are exchangeable, the only thing that matters is the "pattern" of which of the indices \(i_{1},j_{1},i_{2},j_{2}\) are repeated vs which are distinct. For example, there will be \(n\) index configurations where \(i_{1}=j_{1}=i_{2}=j_{2}\) are all equal. All \(n\) of these give same contribution. There are \((n)_{4}=n(n-1)(n-2)(n-3)\) configurations where \(i_{1},j_{1},i_{2},j_{2}\) are all distinct. Knowing which indices are repeated/distinct allows us to evaluate the corresponding term in (64). We use the following formal notion of a pattern to organize this idea of repeated vs. distinct indices.

**Definition 17**: _A **pattern** for \((i_{1},j_{1}),(i_{2},j_{2})\) is a subset of all possible index configurations \((i_{1},j_{1}),(i_{2},j_{2})\in[n]^{4}\) represented by an assignment of each index to the letters \(a,b,c,d\). Each letter \(a,b,c,d\) represents a choice of unique indices from \([n]\)._

For example, the pattern \((i_{1},j_{1}),(i_{2},j_{2})=(a,a),(a,a)\) represents the set of all index configurations where all indices are equal and the pattern \((i_{1},j_{1}),(i_{2},j_{2})=(a,b),(c,d)\) represents the set with all indices unique. The pattern \((i_{1},j_{1}),(i_{2},j_{2})=(a,b),(a,c)\) represents all configurations where \(i_{1}=i_{2}\) and \(j_{1},j_{2}\) are unique and different from \(i_{1}=i_{2}\). For this pattern, there are \((n)_{3}=n(n-1)(n-2)\) configurations by filling in \(a,b,c\) with unique numbers in \([n]\). More generally, for a pattern with \(k\) letters, there are \((n)_{k}\) configurations that fall into that pattern.

Fortunately, when enumerating (64), many patterns have _no_ contribution and can be ignored. We formalize this in the following definition.

**Definition 18**: _We say that the configuration of indices \((i_{1},j_{1}),(i_{2},j_{2})\) is **reducible** if \(\{i_{1},j_{1}\}\cap\{i_{2},j_{2}\}=\varnothing\). Otherwise, the index configuration is called **irreducible**. A pattern is called reducible if all index configuration in that pattern are reducible._

By the independence of the random variables \(f_{1}(G_{i_{1}},G_{j_{1}})\) and \(f_{2}(G_{i_{2}},G_{j_{2}})\), whenever \((i_{1},j_{1}),(i_{2},j_{2})\) is reducible, we see that the corresponding term in (64) completely vanishes! Therefore, to evaluate (64), we have only to understand the contribution of _irreducible configurations_. The irreducible configurations can be organized into _irreducible patterns_. For example, the pattern \((a,b),(c,c)\) is reducible (since formally \(\{a,b\}\cap\{c\}=\varnothing\)) and so any configuration from this pattern has _no_ contribution in the expectation.

There are 11 irreducible patterns. (All these patterns are listed as part of Table 4.) The expected value of the terms for each pattern will give a contribution that is expressed in terms of the \(J_{a,b}\) depending on the details of exactly which indices are repeated. Then by enumerating the number of configurations in each pattern, we can evaluate (64). This strategy is precisely how we evaluate each variance/covariance in this section.
Calculation of \(\mathbf{Var}\left[R^{\ell+1}\right]\) :

First, applying the identity in (17), we get

\[\mathbf{Var}\left[R^{\ell+1}\right] =\left(\frac{2}{n_{\ell}}\right)^{4}\mathbf{Var}\left[\sum_{i,j=1} ^{n_{\ell}}\varphi^{2}(G_{i})\varphi^{2}(\hat{G}_{j})\right]\] \[=\frac{16}{n_{\ell}^{4}}\left(\mathbf{E}\left[\sum_{\begin{subarray} {c}i_{1},j_{1}\\ i_{2},j_{2}\end{subarray}}\varphi^{2}(G_{i_{1}})\varphi^{2}(\hat{G}_{j_{1}}) \varphi^{2}(G_{i_{2}})\varphi^{2}(\hat{G}_{j_{2}})\right]-\mathbf{E}\left[ \sum_{i,j=1}^{n_{\ell}}\varphi^{2}(G_{i})\varphi^{2}(\hat{G}_{j})\right]^{2} \right).\]

\(\mathbf{Var}[R^{\ell+1}]\) follows the form of (64), with \(f_{1}(G_{i},\hat{G}_{i})=f_{2}(G_{i},\hat{G}_{i})=\varphi^{2}(G_{i})\varphi^{2 }(\hat{G}_{j})\). We then evaluate the contribution from each irreducible pattern in Table 4. Combining all these cases and simplifying based on powers of \(\frac{1}{n_{\ell}}\), we arrive at the following expression:

\[\mathbf{Var}\left[R^{\ell+1}\right] =\frac{4}{n_{\ell}}(J_{2,2}+1)+\frac{16}{n_{\ell}^{2}}\left(2J_{4,2}-\frac{5}{2}J_{2,2}+J_{2,2}^{2}+\frac{5}{8}\right)\] \[+\frac{16}{n_{\ell}^{3}}\left(J_{4,4}-2J_{4,2}-2J_{2,2}^{2}+2J_{2,2}-\frac{9}{8}\right).\]

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|} \hline \# & \((i_{1},j_{1})\) & \((i_{2},j_{2})\) & \(\mathbf{E}[f_{1}(G_{i_{1}},\hat{G}_{j_{1}})]\) & \(\mathbf{E}[f_{2}(G_{i_{2}},\hat{G}_{j_{2}})]\) & \(\mathbf{E}[f_{1}(G_{i_{1}},\hat{G}_{j_{1}})f_{2}(G_{i_{2}},\hat{G}_{j_{2}})]\) \\ \hline \hline \((n)_{1}\) & \((a,a)\) & \((a,a)\) & \(J_{2,2}\) & \(J_{2,2}\) & \(J_{4,4}\) \\ \hline \multirow{6}{*}{\((n)_{2}\)} & \((a,b)\) & \((a,b)\) & \multirow{2}{*}{\(\left(\frac{1}{2}\right)^{2}\)} & \multirow{2}{*}{\(\left(\frac{1}{2}\right)^{2}\)} & \multirow{2}{*}{\(\left(\frac{3}{2}\right)^{2}\)} \\ \cline{2-3} \cline{5-6}  & \((a,b)\) & \((b,a)\) & & & \(J_{2,2}\) \\ \cline{2-3} \cline{5-6}  & \((a,a)\) & \((a,b)\) & & & \\ \cline{2-3} \cline{5-6}  & \((a,b)\) & \((a,a)\) & & & \(\frac{1}{2}J_{4,2}\) \\ \cline{2-3} \cline{5-6}  & \((b,a)\) & \((a,a)\) & & & \\ \hline \multirow{6}{*}{\((n)_{3}\)} & \((a,b)\) & \((a,c)\) & \multirow{6}{*}{\(\left(\frac{1}{2}\right)^{2}\)} & \multirow{6}{*}{\(\left(\frac{1}{2}\right)^{2}\)} & \multirow{6}{*}{\(\left(\frac{1}{2}\right)^{2}\)} \\ \cline{2-3} \cline{5-6}  & \((a,b)\) & \((c,b)\) & & & \\ \cline{1-1} \cline{2-3} \cline{5-6}  & \((a,b)\) & \((c,a)\) & & & \\ \cline{1-1} \cline{2-3} \cline{5-6}  & \((a,b)\) & \((b,c)\) & & & \\ \hline \end{tabular}
\end{table}
Table 4: \(\mathbf{Var}[R^{\ell+1}]\) calculated in the form of (64) with \(f_{1}(G_{i},\hat{G}_{j})=f_{2}(G_{i},\hat{G}_{j})=\varphi^{2}(G_{i})\varphi^{2 }(\hat{G}_{j})\). The contribution from all 11 possible _irreducible_ patterns of the indices are shown.
Calculation of \(\mathbf{Var}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right]\) :

Applying identity (16), we can express \(\mathbf{Var}[R^{\ell+1}\sin^{2}(\theta^{\ell+1})]\) as

\[\mathbf{Var}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right]\] \[=\frac{1}{4}\left(\frac{2}{n_{\ell}}\right)^{4}\mathbf{Var}\left[ \sum_{i,j=1}^{n_{\ell}}\left(\varphi(G_{i})\varphi(\hat{G}_{j})-\varphi(G_{j}) \varphi(\hat{G}_{i})\right)^{2}\right]\] \[=\frac{4}{n_{\ell}^{4}}\sum_{\begin{subarray}{c}i_{1},j_{1}\\ i_{2},j_{2}\end{subarray}}\mathbf{Cov}\left(\left(\varphi(G_{i_{1}})\varphi( \hat{G}_{j_{1}})-\varphi(G_{j_{1}})\varphi(\hat{G}_{i_{1}})\right)^{2},\left( \varphi(G_{i_{2}})\varphi(\hat{G}_{j_{2}})-\varphi(G_{j_{2}})\varphi(\hat{G}_ {i_{2}})\right)^{2}\right)\]

Note that we can express \(\mathbf{Var}[R^{\ell+1}\sin^{2}(\theta^{\ell+1})]\) as in (64) by letting \(f_{1}(G_{i},\hat{G}_{j})=f_{2}(G_{i},\hat{G}_{j})=(\varphi(G_{i})\varphi( \hat{G}_{j})-\varphi(G_{j})\varphi(\hat{G}_{i}))^{2}\).We then evaluate the contribution from each irreducible pattern in Table 5. Combining all these cases and simplifying based on powers of \(\frac{1}{n_{\ell}}\), we arrive at the following expression:

\[\mathbf{Var}\left[R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right]\] \[=\frac{8}{n_{\ell}}\left(-8J_{1,1}^{4}+8J_{1,1}^{2}J_{2,2}+4J_{1, 1}^{2}-8J_{1,1}J_{3,1}+J_{2,2}+1\right)\] \[+\frac{2}{n_{\ell}^{2}}\left(80J_{1,1}^{4}-96J_{1,1}^{2}J_{2,2}-4 0J_{1,1}^{2}+96J_{1,1}J_{3,1}+24J_{2,2}^{2}-12J_{2,2}-32J_{3,1}^{2}+5\right)\] \[+\frac{2}{n_{\ell}^{3}}\left(-48J_{1,1}^{4}+64J_{1,1}^{2}J_{2,2}+ 24J_{1,1}^{2}-64J_{1,1}J_{3,1}-24J_{2,2}^{2}+8J_{2,2}+32J_{3,1}^{2}-9\right).\]

\(\mathbf{Var}[R^{\ell+1}\sin^{2}(\theta^{\ell+1})]\) Calculation

\begin{tabular}{|c|c|c|c|c|c|} \hline \# & \((i_{1},j_{1})\) & \((i_{2},j_{2})\) & \(\mathbf{E}[f(G_{i_{1}},\hat{G}_{j_{1}})]\) & \(\mathbf{E}[f(G_{i_{2}},\hat{G}_{j_{2}})]\) & \(\mathbf{E}[f(G_{i_{1}},\hat{G}_{j_{1}})f(G_{i_{2}},\hat{G}_{j_{2}})]\) \\ \hline \hline \multirow{2}{*}{\((n)_{2}\)} & \((a,b)\) & \((a,b)\) & & & \multirow{2}{*}{\(6J_{2,2}^{2}-8J_{3,1}^{2}+\frac{9}{2}\)} \\ \cline{2-6}  & \((a,b)\) & \((b,a)\) & & & \\ \hline \multirow{2}{*}{\((n)_{3}\)} & \((a,b)\) & \((a,c)\) & \multirow{2}{*}{\(\frac{1}{2}-2J_{1,1}^{2}\)} & \(\frac{1}{2}-2J_{1,1}^{2}\) & \multirow{2}{*}{\(4J_{2,2}J_{1,1}^{2}-4J_{3,1}J_{1,1}+\frac{1}{2}J_{2,2}+\frac{3}{4}\)} \\ \cline{2-6}  & \((a,b)\) & \((c,b)\) & & & \\ \cline{1-1} \cline{2-6}  & \((a,b)\) & \((b,c)\) & & & \\ \hline \end{tabular} Table 5: \(\mathbf{Var}[R^{\ell+1}\sin^{2}(\theta^{\ell+1})]\) calculated in the form of (64) with \(f_{1}(G_{i},\hat{G}_{j})=f_{2}(G_{i},\hat{G}_{j})\)\(=(\varphi(G_{i})\varphi(\hat{G}_{j})-\varphi(G_{j})\varphi(\hat{G}_{i}))^{2}\). The _non-zero_ contribution _irreducible_ patterns of the indices are shown. Note that because \(f_{1}(G_{i_{1}},G_{j_{1}})=0\) when \(i_{1}=j_{1}\) and \(f_{2}(G_{i_{2}},G_{j_{2}})=0\) when \(i_{2}=j_{2}\), there are 5 irreducible patterns (of the possible 11) that have zero contribution and are not displayed in this table.
\[\mathbf{Cov}\left(R^{\ell+1}\sin^{2}(\theta^{\ell+1}),R^{\ell+1}\right)=\mathbf{E} \left[\left(R^{\ell+1}\right)^{2}\sin^{2}(\theta^{\ell+1})\right]-\mathbf{E} \left[R^{\ell+1}\sin^{2}(\theta^{\ell+1})\right]\mathbf{E}\left[R^{\ell+1} \right].\]

Applying known identities (16, 17) derived in Appendix A.5 and Appendix A.6, we can express this in the form of (64), where \(f_{1}(G_{i},\hat{G}_{j})=(\varphi(G_{i})\varphi(\hat{G}_{j})-\varphi(G_{j}) \varphi(\hat{G}_{i}))^{2}\), and \(f_{2}(G_{i},\hat{G}_{j})=\varphi^{2}(G_{i})\varphi^{2}(\hat{G}_{j})\). Table 6 shows the calculation of all the irreducible patterns. Collecting all cases and simplifying based on powers of \(\frac{1}{n_{\ell}}\) gives:

\[\mathbf{Cov}\left(R^{\ell+1}\sin^{2}(\theta^{\ell+1}),R^{\ell+1} \right)=\frac{1}{n_{\ell}}\left(16J_{1,1}^{2}-32J_{1,1}J_{3,1}+8J_{2,2}+8\right)\] \[+\frac{1}{n_{\ell}^{2}}\left(32J_{1,1}^{2}J_{2,2}-40J_{1,1}^{2}+ 96J_{1,1}J_{3,1}-32J_{1,1}J_{3,3}+16J_{2,2}^{2}-32J_{2,2}-32J_{3,1}^{2}+16J_{ 4,2}+10\right)\] \[+\frac{1}{n_{\ell}^{3}}\left(24J_{1,1}^{2}-32J_{1,1}^{2}J_{2,2}-64 J_{1,1}J_{3,1}+32J_{1,1}J_{3,3}-16J_{2,2}^{2}+24J_{2,2}+32J_{3,1}^{2}-16J_{ 4,2}-18\right).\]

### Infinite Width Update Rule

**Lemma 19**: _Let \(f(x)\) be a feed forward neural network as defined in 2, and assume each layer has infinite width (ie. \(n_{\ell}=\infty\) for \(1\leq\ell\leq L\)). Then, the angle \(\theta^{\ell}\) between inputs at layer \(\ell\) of \(f\) follows the update rule given by_

\[\cos(\theta^{\ell+1})=2J_{1,1}(\theta^{\ell}).\]
**Proof** We begin with the identity in (15), and use the inner product to introduce a \(\cos(\theta^{\ell+1})\) term into our equation.

\[\frac{\|\varphi_{\alpha}^{\ell}\|\|\varphi_{\beta}^{\ell}\|}{n_{\ell}}\sum_{i=1} ^{n_{\ell}}2\varphi(G_{i})\varphi(\hat{G}_{i})=\langle\varphi_{\alpha}^{\ell+1},\varphi_{\beta}^{\ell+1}\rangle=\|\varphi_{\alpha}^{\ell+1}\|\|\varphi_{ \beta}^{\ell+1}\|\cos(\theta^{\ell+1}).\]

Applying the identities in (14) to \(\|\varphi_{\alpha}^{\ell+1}\|\) and \(\|\varphi_{\beta}^{\ell+1}\|\), we get

\[\frac{\|\varphi_{\alpha}^{\ell}\|\|\varphi_{\beta}^{\ell}\|}{n_{ \ell}}\sum_{i=1}^{n_{\ell}}2\varphi(G_{i})\varphi(\hat{G}_{i}) =\sqrt{\frac{\|\varphi_{\alpha}^{\ell+1}\|^{2}}{n_{\ell}}\sum_{i= 1}^{n_{\ell}}2\varphi^{2}(G_{i})}\sqrt{\frac{\|\varphi_{\beta}^{\ell+1}\|^{2} }{n_{\ell}}\sum_{i=1}^{n_{\ell}}2\varphi^{2}(\hat{G}_{i})}\cos(\theta^{\ell+1})\] \[\frac{1}{n_{\ell}}\sum_{i=1}^{n_{\ell}}\varphi(G_{i})\varphi(\hat {G}_{i}) =\sqrt{\frac{1}{n_{\ell}}\sum_{i=1}^{n_{\ell}}\varphi^{2}(G_{i}) }\sqrt{\frac{1}{n_{\ell}}\sum_{i=1}^{n_{\ell}}\varphi^{2}(\hat{G}_{i})}\cos( \theta^{\ell+1}).\]

Taking the limit as \(n_{\ell}\to\infty\), and using the Law of Large Numbers, we have

\[\lim_{n_{\ell}\to\infty}\left(\frac{1}{n_{\ell}}\sum_{i=1}^{n_{ \ell}}\varphi(G_{i})\varphi(\hat{G}_{i})\right) =\lim_{n_{\ell}\to\infty}\left(\sqrt{\frac{1}{n_{\ell}}\sum_{i=1} ^{n_{\ell}}\varphi^{2}(G_{i})}\sqrt{\frac{1}{n_{\ell}}\sum_{i=1}^{n_{\ell}} \varphi^{2}(\hat{G}_{i})}\cos(\theta^{\ell+1})\right)\] \[\mathbf{E}\left[\varphi(G_{i})\varphi(\hat{G}_{i})\right] =\sqrt{\mathbf{E}\left[\varphi^{2}(G_{i})\right]}\sqrt{\mathbf{E }[\varphi^{2}(\hat{G}_{i})]}\cos(\theta^{\ell+1})\] \[J_{1,1}(\theta^{\ell}) =\frac{1}{2}\cos(\theta^{\ell+1})\implies\cos(\theta^{\ell+1})=2J _{1,1}(\theta^{\ell}).\]

Note that in the \(J_{p}\) notation from Cho and Saul (2009), the translation is \(2J_{1,1}=\frac{J_{1}}{\pi}\), so this is is the result of the theorem.

**Appendix B.**

### Derivation of Lower-Order J Functions - Proof of Proposition 2

**Proof** [Of formula for \(J_{0,0}\)] We find a differential equation that \(J_{0,0}\) satisfies and solve it to obtain the formula. First note the initial condition \(J_{0,0}(0)=\mathbf{E}[1\{G>0\}]=\frac{1}{2}\). To find \(J_{0,0}^{\prime}(\theta)\), we take the derivative into the expectation and have by the chain rule that

\[J_{0,0}^{\prime}(\theta) =\mathbf{E}[1\{G>0\}1^{\prime}\{G\cos\theta+W\sin\theta>0\}G](- \sin\theta)\] \[+\mathbf{E}[1\{G>0\}1^{\prime}\{G\cos\theta+W\sin\theta>0\}W]\cos\theta.\]

Applying the change of variables as in (37), we have

\[J_{0,0}^{\prime}(\theta) =\mathbf{E}[(Z\cos\theta+Y\sin\theta)1\{Z\cos\theta+Y\sin\theta>0 \}1^{\prime}\{Z>0\}](-\sin\theta)\] \[+\mathbf{E}[(Z\sin\theta-Y\cos\theta)1\{Z\cos\theta+Y\sin\theta>0 \}1^{\prime}\{Z>0\}]\cos\theta\] \[=\mathbf{E}[(Y\sin\theta)1\{Y\sin\theta>0\}]\frac{-\sin\theta}{ \sqrt{2\pi}}+\mathbf{E}[(-Y\cos\theta)1\{Y\sin\theta>0\}]\frac{\cos\theta}{ \sqrt{2\pi}}\] \[=(-\sin^{2}\theta-\cos^{2}\theta)\mathbf{E}[Y1\{Y>0\}]\frac{1}{ \sqrt{2\pi}}=-\frac{1}{2\pi},\]
where we have used (32) to evaluate the integrals involving \(1^{\prime}\{Z>0\}\) and \(\mathbf{E}[Y1\{Y>0\}]=(\sqrt{2\pi})^{-1}\) from Lemma 7. We now have \(J^{\prime}_{0,0}(\theta)=-\frac{1}{2\pi}\) with initial condition given by \(J_{0,0}(0)=0\). Solving this differential equation gives the desired result. \(\blacksquare\)

**Proof** [\(J_{1,0}\) and \(J_{1,1}\)] Here we use the Gaussian integration-by-parts strategy (31 -32).

**Formula for \(J_{1,0}\):**

\[J_{1,0}(\theta) =\mathbf{E}[G1\{G>0\}\ 1\{G\cos\theta+W\sin\theta>0\}]\] \[=\mathbf{E}\left[\frac{d}{dg}\left(\ 1\{G>0\}\ 1\{G\cos\theta+W\sin \theta>0\}\right)\right]\] \[=\mathbf{E}\left[1^{\prime}\{G>0\}\ 1\{G\cos\theta+W\sin\theta>0\} \right]+\mathbf{E}\left[1\{G>0\}\ 1^{\prime}\{G\cos\theta+W\sin\theta>0\}\right]\cos\theta.\]

By using the change of variables as in (37) on the second term, we arrive at

\[J_{1,0}(\theta) =\mathbf{E}[1\{W\sin\theta>0\}]\frac{1}{\sqrt{2\pi}}+\cos\theta \mathbf{E}[1\{Z\cos\theta+Y\sin\theta>0\}1^{\prime}\{Z>0\}]\] \[=\frac{1}{2}\frac{1}{\sqrt{2\pi}}+\cos\theta\mathbf{E}[1\{Y\sin \theta>0\}]\frac{1}{\sqrt{2\pi}}=\frac{1}{2}\frac{1}{\sqrt{2\pi}}+\frac{\cos \theta}{2}\frac{1}{\sqrt{2\pi}}=\frac{1+\cos\theta}{2\sqrt{2\pi}}.\]

**Formula for \(J_{1,1}\):**

\[J_{1,1}(\theta) =\mathbf{E}[G(G\cos\theta+W\sin\theta)1\{G>0\}1\{G\cos\theta+W\sin \theta>0\}]\] \[=\mathbf{E}\left[\frac{d}{dg}((G\cos\theta+W\sin\theta)1\{G>0\}1 \{G\cos\theta+W\sin\theta>0\})\right]\] \[=\mathbf{E}[\cos\theta\ 1\{G>0\}1\{G\cos\theta+W\sin\theta\}]\] \[+\mathbf{E}[(G\cos\theta+W\sin\theta)1^{\prime}\{G>0\}1\{G\cos \theta+W\sin\theta>0\}]\] \[+\mathbf{E}[(G\cos\theta+W\sin\theta)1\{G>0\}1^{\prime}\{G\cos \theta+W\sin\theta>0\}]\cos\theta\] \[=\cos\theta J_{0,0}+\mathbf{E}[W\sin\theta\ 1\{W\sin\theta>0\}] \frac{1}{\sqrt{2\pi}}+\mathbf{E}[Z1\{Z\cos\theta+Y\sin\theta\}1^{ \prime}\{z>0\}]\] \[=\cos\theta J_{0,0}+\sin\theta\mathbf{E}[\varphi(W)]\frac{1}{ \sqrt{2\pi}}+0=\frac{\sin\theta+(\pi-\theta)\cos\theta}{2\pi}.\]

\(\blacksquare\)

### Proof of Explicit Formulas for \(J_{n,0}\) and \(J_{n,1}\)

Once the recursion is established, the formula for both \(J_{n,0}\) and \(J_{n,1}\) is a simple proof by induction. We provide a detailed proof for \(J_{n,0}\) here; \(J_{n,1}\) is similar.

**Lemma 20**: _Let \(J_{n,0}^{rec}\) be the recursively defined formula, and let \(J_{n,0}^{exp}\) is the explicitly defined formula for \(J_{n,0}\), namely_

\[J_{n,0}^{rec}:=(n-1)J_{n-2,0}^{rec}+\frac{\sin^{n-1}\theta\cos\theta}{c_{n\bmod 2 }}(n-2)!!\,,\quad J_{1,0}^{rec}:=J_{1,0},\quad J_{0,0}^{rec}:=J_{0,0}\]
\[J_{n,0}^{exp}:=(n-1)!!\left(J_{n\bmod 2,0}+\frac{\cos\theta}{c_{n\bmod 2}}\sum_{ \begin{subarray}{c}i\neq n(\bmod 2)\\ 0<i<n\end{subarray}}\frac{(i-1)!!}{i!\,!}\sin^{i}\theta\right).\]

_Then \(J_{n,0}^{rec}=J_{n,0}^{exp}\) for all \(n\geq 0\)._

**Proof** Let \(S_{n}\), \(n\in\mathbb{N},\;n\geq 2\) be the statement \(J_{n,0}^{rec}=J_{n,0}^{exp}\) and \(J_{n-1,0}^{rec}=J_{n-1,0}^{exp}\). We prove \(S_{n}\) is true by induction. The base case \(S_{2}\) is true because,

\[J_{2,0}^{rec} =(2-1)J_{0,0}+\frac{\sin\theta\cos\theta}{c_{2\bmod 2}}(2-2)!!=J_{ 0,0}+\frac{\cos\theta\sin\theta}{2\pi}\] \[J_{2,0}^{exp} =(2-1)!!\,J_{0,0}+\cos\theta\sum_{i=1}^{1}\frac{(2-1)!!}{(2i-1)!!}(2i-2)!!\,\frac{\sin^{2i-1}\theta}{2\pi}=J_{0,0}+\frac{\cos\theta\sin\theta }{2\pi},\]

and the fact that \(J_{1,0}^{rec}=J_{1,0}^{exp}\) is trivial. Induction step: Assume \(S_{n}\) is true. To prove \(S_{n+1}\), we have only to show that \(J_{n+1,0}^{exp}=J_{n+1,0}^{rec}\). To do this, we separate the last term of the sum to get

\[J_{n+1,0}^{exp}= n!\,!\left(J_{(n+1)\bmod 2,0}+\frac{\cos\theta}{c_{(n+1)\bmod 2} }\sum_{\begin{subarray}{c}i\neq(n+1)(\bmod 2)\\ 0<i<n-1\end{subarray}}\frac{(i-1)!!}{i!\,!}\sin^{i}\theta\right)\] \[+n!\,!\,\frac{\cos\theta}{c_{(n+1)\bmod 2}}\frac{(n-1)!!}{n!\,!} \sin^{n}\theta.\]

Because the parity of \(n+1\) and \(n-1\) are the same, and using \(n!\,!=n(n-2)!\,!\) we recognize the first term as \(nJ_{n-1,0}^{exp}\). So after simplifying the last term, we remain with

\[J_{n+1,0}^{exp}=nJ_{n-1,0}^{exp}+\frac{\sin^{n}\theta\cos\theta}{c_{(n+1) \bmod 2}}(n-1)!!=J_{n+1,0}^{rec},\]

by the induction hypothesis. This completes the induction.

### Bijection between Paths in Graphs of J Functions and the Bessel Number graphs \(P\),\(Q\)

Let \(G_{J^{*}}=(V_{J*},E_{J*})\) be the graph of \(J_{a,b}^{*}\) as in Figure 5c. Similarly, let \(G_{P}=(V_{P},E_{P})\) and \(G_{Q}=(V_{Q},E_{Q})\) be the graph of the \(P\) and \(Q\) matrices up to row \(a\), respectively, as in Figures 5a, 5b. We define a map \(\lambda:\mathbb{Z}^{2}\times\mathbb{Z}^{2}\to\mathbb{Z}^{2}\) as follows: Let \(((i,j),(m,n))\in\mathbb{Z}^{2}\times\mathbb{Z}^{2},\;0\leq i\leq a,\;b-a+m\leq j \leq b\). Then define \(\lambda\) by

\[\lambda((i,j),(m,n)):=(i-m,j-n),\quad\lambda^{-1}((i,j),(m,n))=(i+m,j+n).\]

The function \(\lambda\) can be used as a map between vertices of graph \(G_{J^{*}}\) to vertices of graph \(G_{P}\) or \(G_{Q}\). Let \(\pi=(v_{1},v_{2},...,v_{k-1},v_{k})\) be a path in \(G_{J^{*}}\) from vertex \(v_{1}=(m,n)\) to vertex
\(v_{k}=(a,b)\), where \(v_{i}\in\mathbb{Z}^{2},\ 1\leq i\leq k\) is a vertex on the graph. \(\lambda\) extends to a map on paths, \(\Lambda\), defined by

\[\Lambda((v_{1},v_{2},...,v_{k-1},v_{k})):= (\lambda(v_{1},v_{1}),\lambda(v_{2},v_{1}),...,\lambda(v_{k-1},v_{1 }),\lambda(v_{k},v_{1})),\] \[\Lambda^{-1}((v_{1},v_{2},...,v_{k-1},v_{k}))= (\lambda^{-1}(v_{1},v_{1}),\lambda^{-1}(v_{2},v_{1}),...,\lambda^{- 1}(v_{k-1},v_{1}),\lambda^{-1}(v_{k},v_{1})).\]

Now, let \(\Gamma_{J^{*}}(a,b,m,n)\) be the set of all paths in the graph of \(J^{*}\) from \(J^{*}_{m,n}\) to \(J^{*}_{a,b}\), and let \(\Gamma_{P}(a,b,m,n)\) be the set of all paths in the graph of \(P\)\(P(0,0)\) to \(P(a-m,b-n)=P(\lambda((a,b),(m,n)))\). For example, \(\Gamma_{J^{*}}(6,8,0,4)\) is the set of all paths which run from \(J^{*}_{6,8}\) to \(J^{*}_{0,4}\), and \(\Gamma_{P}(6,8,0,4)\) is the set of all paths which run from \(P(0,0)\) to \(P(6,4)\).

With these definitions, \(\Lambda:\Gamma_{J^{*}}(a,b,0,n)\to\Gamma_{P}(a,b,0,n)\) is a bijection. An illustration of all paths \(\pi\in\Pi(6,8,0,6)\) and the corresponding paths \(\Lambda(\pi)\in\Gamma_{P}(6,8,0,6)\) is given in Figure 6. Similarly, if we let \(\Gamma_{Q}(a,b,m,n)\) be the set of all paths from \(Q(0,0)\) to \(Q(a-m,b-n)=Q(\lambda((a,b),(m,n)))\) then \(\Lambda:\Gamma_{J^{*}}(a,b,1,n)\to\Gamma_{Q}(a,b,1,n)\) is a bijection. This bijection establishes the equality of the weighted paths claim in (48).

Figure 6: Top: All paths \(\pi\in\Gamma_{J^{*}}(6,8,0,6)\). Bottom: All paths \(\Lambda(\pi)\in\Gamma_{P}(6,8,0,6)\). The paths are lined up so that for each path \(\pi\) in the top row, \(\Lambda(\pi)\) appears in the bottom row.
## Appendix C Recursions for the \(P\) and \(Q\) numbers - Proof of Lemma 10

Earlier work established the following properties of the \(P\) and \(Q\) numbers

**Theorem 21** (Kreinin (2016)): _The elements of the matrices \(P\) and \(Q\) satisfy_

\[b\cdot P(a,b) =a\cdot P(a-1,b-1), a\geq 1,1\leq b\leq a \tag{66}\] \[P(a+1,b) =P(a,b-1)+(b+1)\cdot P(a,b+1), a\geq 0,1\leq b\leq a\] (67) \[Q(a,b) =P(a,b)+(b+1)\cdot Q(a-1,b+1), a\geq 1,1\leq b\leq a\] (68) \[Q(a,b) =a\cdot Q(a-2,b)+Q(a-1,b-1), a\geq 2,1\leq b\leq a. \tag{69}\]

[Of Lemma 10] Equation (67) tells us that \(P(a,b)=P(a-1,b-1)+(b+1)\cdot P(a-1,b+1)\) for \(a\geq 1,1\leq b\leq a-1\), while equation (66) tells us that \(P(a-1,b+1)=\frac{(a-1)}{(b+1)}\cdot P(a-2,b)\) for \(a\geq 2,0\leq b\leq a-2\). Putting these together, we get the following recurrence equation for \(P(a,b)\):

\[P(a,b) =P(a-1,b-1)+(b+1)\left(\frac{(a-1)}{(b+1)}\cdot P(a-2,b)\right)\] \[=(a-1)\cdot P(a-2,b)+P(a-1,b-1),\]

which holds for \(a\geq 3,1\leq b\leq a-2\). Further, looking at equation (69), we see that the recursion for the \(Q\) numbers is very similar to that of the \(P\) numbers, but with a coefficient of \(a\) rather than \((a-1)\). This establishes Lemma 10.
