Breaking Neural Network Scaling Laws with Modularity

## Abstract

Abstract Modular neural networks outperform nonmodular neural networks on tasks ranging from visual question answering to robotics. These performance improvements are thought to be due to modular networks’ superior ability to model the compositional and combinatorial structure of real-world problems. However, a theoretical explanation of how modularity improves generalizability, and how to leverage task modularity while training networks remains elusive. Using recent theoretical progress in explaining neural network generalization, we investigate how the amount of training data required to generalize on a task varies with the intrinsic dimensionality of a task’s input. We show theoretically that when applied to modularly structured tasks, while nonmodular networks require an exponential number of samples with task dimensionality, modular networks’ sample complexity is independent of task dimensionality: modular networks can generalize in high dimensions. We then develop a novel learning rule for modular networks to exploit this advantage and empirically show the improved generalization of the rule, both in- and out-of-distribution, on high-dimensional, modular tasks.

## 1 Introduction

Modular neural network (NN) architectures have achieved impressive results in a variety of domains ranging from visual question answering (VQA) (Andreas et al., ; ; Hu et al., ; Johnson et al., ; Yi et al., ; Kim et al., ), reinforcement learning (Goyal et al., ; Madan et al., ), robotics (Alet et al., ; Pathak et al., ; Yang et al., ) and natural language processing for which modular architectures based on attention (Bahdanau et al., ) are standard. Modular NNs are thought to be better generalized by facilitating combinatorial generalization, a phenomenon where a learning system recombines previously learned components in novel ways to generalize to unseen task inputs (Alet et al., ; D’Amario et al., ; Mittal et al., ; ; Jarvis et al., ). Yet, a fundamental understanding of why modularity benefits generalization is lacking.

In parallel, the generalization properties of monolithic (nonmodular) NNs have been increasingly well understood both theoretically and empirically. In particular, current theory can explain the double descent phenomenon where NN generalization error decreases with increasingly large capacity (Belkin et al., ; Spigler et al., ; Neal et al., ; Rocks & Mehta, ).
NN learning in a certain regime can also be understood as kernel regression (Jacot et al., ). Extensive empirical studies have also measured scaling laws of NN generalization error (Kaplan et al., ), and moreover, these scaling laws can be explained theoretically (Bahri et al., ; Hutter, ; Hastie et al., ). However, these laws indicate that the sample complexity required to generalize on a task scales exponentially with the intrinsic dimensionality of the task’s input (McRae et al., ; Sharma & Kaplan, ). This raises the question: how can we hope to generalize on high-dimensional problems with limited training data?

In this work, we investigate how modular NNs can circumvent this exponential number of samples. We first synthesize existing generalization results in a simple theoretical model of NN generalization error and empirically validate it on tasks with varying intrinsic dimensionality. We then use our model to show theoretically that appropriately structured modular NNs avoid using an exponential number of samples on modular tasks. However, recent work shows that architectural modularity is in practice not sufficient on its own to solve modular tasks efficiently (Csordás et al., ; Mittal et al., ); a solution to align NN modules to a task’s modularity is lacking. We propose a novel learning rule that aligns NN modules to approach the underlying modules of a task and empirically demonstrate its improved generalization.

We summarize our contributions as follows:

- •
We propose a simple, theoretical model of NN learning that synthesizes existing generalization results. Our model predicts the generalization error of NNs under varying number of model parameters, number of training samples, and dimensions of variation in a task input. We empirically validate our theoretical model on a novel parametrically controllable sine wave regression task and show that sample complexity varies exponentially with task dimension.
- •
We apply the theoretical model to compute explicit, non-asymptotic expressions for generalization error in modular architectures; to our knowledge, we are the first to do so. Our result demonstrates that sample complexity is independent of task dimension for modular NNs applied to modular tasks of a specific form.
- •
Based on our theory, we develop a learning rule to align NN modules to the modules underlying high-dimensional modular tasks with the goal of promoting generalization on these tasks.
- •
We empirically validate the improved generalizability (both in- and out-of-distribution) of our modular learning approach on parametrically controllable, high-dimensional tasks: sine-wave regression and Compositional CIFAR-10.

## 2 Related Work

### 2.1 Modular Neural Networks

Recent efforts to model cognitive processes show that functional modules and compositional representations emerge after training on a task (Yang et al., ; Yamashita & Tani, ; Iyer et al., ). Partly inspired by this, recent works in AI propose using modular networks: networks composed of sparsely connected, reusable modules (Alet et al., ; ; Chang et al., ; Chaudhry et al., ; Shazeer et al., ; Ashok et al., ; Yang et al., ; Sax et al., ; Pfeiffer et al., ). Empirically, modularity improves out-of-distribution generalization (Bengio et al., ; Madan et al., ; Mittal et al., ; Jarvis et al., ), modular generative models are effective unsupervised learners (Parascandolo et al., ; Locatello et al., ) and modular architectures can be more interpretable (Agarwal et al., ). In addition, meta-learning algorithms can discover and learn the modules without prespecifying them (Chen et al., ; Sikka et al., ; Chitnis et al., ).

Recent empirical studies have investigated how modularity influences network performance and generalization. The degree of modularity increases systematic generalization performance in VQA tasks (D’Amario et al., ) and sequence-based tasks (Mittal et al., ). Rosenbaum et al. () and Cui & Jaech () study routing networks (Rosenbaum et al., ), a type of modular architecture, and identified several difficulties with training these architectures including training instability and module collapse. Csordás et al. () and Mittal et al. () extend this type of analysis to more general networks to show that NN modules often may not be optimally used to promote task performance despite having the potential to do so. These analyses are primarily empirical; in contrast, in our work, we aim to provide a theoretical basis for how modularity may improve generalization. Moreover, given that architectural modularity may not be sufficient to ensure generalization, we propose a learning rule designed to align NN modules to the modularity of the task.

### 2.2 Neural Network Scaling Laws

Many works present frameworks to quantify scaling laws that map a NN’s parameter count or training dataset size to an estimated testing loss. Empirically and theoretically, these works find that testing loss scales as a power-law with respect to the dataset size and parameter count on well-trained NNs (Bahri et al., ; Rosenfeld et al., ), including transformer-based language models (Sharma & Kaplan, ; Clark et al., ; Tay et al., ).

Many previous works also conclude that generalizations of power-law or nonpower-law-based distributions can also model neural scaling laws well, in many cases better than vanilla power-law frameworks (Mahmood et al., ; Alabdulmohsin et al., ). For instance, Hutter () shows that countably infinite parameter models closely follow non-power-law-based distributions under unbounded data complexity regimes. In another case, Sorscher et al. () show that exponential scaling works better than power-law scaling if the testing loss is associated with a pruned dataset size, given a pruning metric that discards easy or hard examples under abundant or scarce data guarantees, respectively.

Some works approach this problem by modeling NN learning as manifold or kernel regression. For example, McRae et al. () considers regression on manifolds and concludes that sample complexity scales based on the intrinsic manifold dimension of the data. In another case, Canatar et al. () draws correlations between the study of kernel regression to how infinite-width deep networks can generalize based on the size of the training dataset and the suitability of a particular kernel for a task. Along these lines, several works use random matrix theory to derive scaling laws for kernel regression (Hastie et al., ; Cui et al., ; ; Wei et al., ; Jin et al., ).

Among other observations, this body of work shows that in the absence of strong inductive biases, high-dimensional tasks have sample complexity growing roughly exponentially with the intrinsic dimensionality of the data manifold. In this work, we borrow the theoretical techniques from this line of work to investigate if learning the modular structure of modular tasks will reduce the sample complexity of training.

## 3 Modeling Neural Network Generalization

**Table 1: Table of symbols.**
| Symbol | Meaning |
| --- | --- |
| $x$ | Task input |
| $\varphi(x)$ | Feature matrix |
| $y(x)$ | Desired task output |
| $W$ | Weights of target function |
| $\hat{y}(x)$ | Output of model |
| $\theta$ | Weights of model |
| $\Lambda$ | Covariance matrix of $W$ |
| $\lambda_{i}$ | Element of $\Lambda$ |
| $n$ | # of training samples |
| $p$ | # of model parameters |
| $P$ | # of total features |
| $d$ | Task output dimensionality |
| $m$ | Task intrinsic dimensionality |
| $u_{i}$ | Module projection vector |
| $U_{i}$ | Module projection matrix |

In this section, we present a toy model of NN learning that treats NNs as linear functions of their parameters; this is along with the lines of prior work such as Bahri et al. (); Canatar et al. (). Although this common theoretical assumption does not directly apply to practical, non-linear architectures, the assumption provides analytical tractability and, moreover, can be shown to predict generalization even in nonlinear networks (e.g. in the Neural Tangent Kernel literature (Jacot et al., )). Our specific analytical approach follows that of Hastie et al. (). Under this setting, we find exact closed-form expressions for expected training and test loss, representing a simplified version of the results in Hastie et al. (). We find that our toy model captures key features of NN generalization applied to a sine wave regression task. We summarize our notation in Tab 1.

### 3.1 Model Setup

We defer the full details of our theoretical model setup to App A. Here we present an overview: we consider a regression task with input $x\in\mathbb{R}^{m}$ and a feature matrix $\varphi(x)\in\mathbb{R}^{d\times P}$, where $m$ is input dimensionality, $d$ is output dimensionality and $P$ is the number of features. We assume that over the data distribution, the features are distributed i.i.d. from a unit Gaussian: $\varphi(x)_{i,j}\sim\mathcal{N}(0,1)$. We consider the limit when $P\to\infty$. Suppose that our regression target function $y:\mathbb{R}^{m}\to\mathbb{R}^{d}$ is constructed linearly from $\varphi(x)$: $y(x)=\varphi(x)W$. Assume $\mathbb{E}[W]=0$ and $\mathbb{E}[WW^{T}]=\Lambda$, where $\Lambda$ is diagonal and $\operatorname{Tr}(\Lambda)$ is finite. Suppose we aim to approximate $y(x)$ using a model $\hat{y}(x)$ constructed as follows: $\hat{y}(x)=\varphi(x)\begin{bmatrix}I\\
0\end{bmatrix}\theta$,
where $\theta\in\mathbb{R}^{p\times 1}$ are model parameters and $p$ is the number of parameters. This corresponds to the model only being able to control $p$ of the $P$ true underlying parameters in the construction of $y$. We decompose $\Lambda$ blockwise as: $\Lambda=\begin{bmatrix}\Lambda_{1}&0\\
0&\Lambda_{2}\end{bmatrix}$
where $\Lambda_{1}\in\mathbb{R}^{p\times p}$ and $\Lambda_{2}\in\mathbb{R}^{(P-p)\times(P-p)}$. To capture the dependence of the target function on the input dimensionality $m$, we parameterize $\Lambda$ as having individual elements $\lambda_{i}=c\left[i^{-\Omega^{-m}}-(i+1)^{-\Omega^{-m}}\right]$. This corresponds to the number of effective dimensions of variation of $W$ scaling exponentially with $m$, which is consistent with prior work (McRae et al., ). We also produce versions of our theoretical results without setting a specific form for $\lambda_{i}$. Finally, we consider learning $\theta$ as the minimum norm interpolating solution.

Figure: Figure 1: Expected training (left) and test (right) set error in a toy model of NN generalization as a function of the number of samples $n$ and the number of model parameters $p$. The output dimensionality is set as $d=1$.
Refer to caption: extracted/5841346/figs/train_error_heatmap_annotated.png

### 3.2 Theoretical Properties

Next, we theoretically analyze the expected training and test set error of the above model.

###### Theorem 1 .

Given a target function $y$ and model $\hat{y}$ estimated as described above, in the limit that $P\to\infty$, the expected test loss when averaging over $x$ and $W$ is:

$$ $\lim_{P\to\infty}\mathbb{E}\left[|y(x)-\hat{y}(x)|^{2}\right]\\ =d\operatorname{Tr}(\Lambda_{2})F(dn,p)-d\frac{\min(dn,p)}{p}\operatorname{Tr} (\Lambda_{1})+d\operatorname{Tr}(\Lambda)\\ =dF(dn,p)c(p+1)^{-\Omega^{-m}}-d\frac{\min(dn,p)}{p}(c-c(p+1)^{-\Omega^{-m}})+dc$ (1) $$

The expected training loss is:

$$ $\lim_{P\to\infty}\frac{1}{n}\mathbb{E}\left[\left\|y(X)-\hat{y}(X)\right\|_{2} ^{2}\right]\\ =\operatorname{Tr}(\Lambda_{2})(dn-\min(dn,p))\\ =\frac{dn-\min(dn,p)}{n}c(p+1)^{-\Omega^{-m}}$ (2) $$

with $F(n,p)$ defined as $F(n,p)=\mathbb{E}\left[\left\|R^{\dagger}\right\|_{F}^{2}\right]$
where $R\in\mathbb{R}^{n\times p}$ has elements drawn i.i.d. from $\mathcal{N}(0,1)$.

Please see App B for a proof and App D for more details on $F(n,p)$. Under general $\lambda_{i}$, training and test error grow with $Tr(\Lambda_{2})$; the specific rate at which they grow or shrink with parameters $m$ and $p$ depends on how rapidly $Tr(\Lambda_{2})$ grows with $m$ and shrinks with $p$. Under the specific parameterization for $\lambda_{i}$ in Eqn 19, Fig 1 plots the value of the training and test set error for varying $n$ and number of parameters, holding $d=1$. Observe that there is a clear interpolation threshold at $p=n$ where the training loss becomes zero; for $p\geq n$ the model has sufficient capacity to perfectly interpolate the training set. The training loss is positive and increases with $n$ in the underparameterized regime ($p<n$) since the model lacks the capacity to fit increasing amounts of data $n$. At the $p=n$ threshold, the test loss dramatically increases, then decreases as $p$ increases beyond $n$. This is consistent with empirically observed behavior of overparameterized NNs (Belkin et al., ). Similarly, test loss decreases as $n$ increases beyond $p$.

Figure: Figure 2: Empirical trends of training (blue) and test (orange) loss over four parametric variations for a NN trained on a sine wave regression task. The parameters varied are: $k$ (number of modules), $m$ (input dimensionality), $p$ (model size) and $n$ (training set size). In the first two plots, each line indicates a different model architecture, and in the last two plots, each line indicates a different choice of $m$ between $5$ and $9$, with $n/p$ fixed at $1000/1153$ respectively (left/right). The light lines are averaged over all other parameters, and bold lines show averages over the light lines. Dashed lines show theoretical predictions.
Refer to caption: extracted/5841346/figs/empirical_sample_complexity_k.png

###### Theorem 1 .

### 3.3 Empirical Validation

Figure: Figure 3: Theoretically predicted trend of $m$ vs. $n$ to achieve a test loss of $1.2$ on a sine wave regression task. Each line indicates a different fully connected NN with a different width and depth. $m$ increases approximately exponentially with $n$.
Refer to caption: extracted/5841346/figs/theoretical_m_vs_n.png

We empirically validate our theoretical model on a parametrically-variable modular sine wave regression task with targets constructed as $y(x)=\frac{1}{\sqrt{k}}\sum_{i=1}^{k}\sum_{j=1}^{\tau}a_{ij}\sin(\omega_{ij}u_
{i}^{T}x+\phi_{ij})$, where $x\in\mathbb{R}^{m}$ are inputs, $y(x)\in\mathbb{R}$ are outputs, $a_{ij},\omega_{ij},\phi_{ij}\in\mathbb{R}$, $u_{i}\in\mathbb{R}^{m}$ are parameters chosen randomly for each target function and $k$ and $\tau$ are fixed (see App E for further task details). We train fully connected ReLU-activated NNs of varying depth and width on the task. The task allows us to quantify how NN generalization depends on a number of factors such as the dimensionality $m$ of the task input, the number of model parameters $p$, the number of samples $n$ and the number of modules $k$ in the construction of the target function. In Fig 2, we find that our theoretical model matches empirical trends of neural network training and test error (see App A for full details). Nevertheless, we note two key discrepancies between empirical and predicted trends: first, the test loss is empirically larger than predicted under low training data. We hypothesize that this may be because of difficulty optimizing for small $n$: indeed, we find that the training loss is larger than expected for small $n$ (in the overparameterized regime ($n<p$), we expect a training loss of 0 0). Second, the error spike at the interpolation threshold is smaller than theoretically predicted. This again may be due to incomplete optimization, given that the interpolation threshold spike can be viewed as highly adverse fitting to spurious training set patterns.

Our theoretical model predicts that the number of training samples required for generalization to a fixed error rate scales exponentially with task dimensionality $m$ (Fig 3). This raises a practical challenge for high-dimensional problems, which may require massive amounts of training data. Next, we aim to circumvent this exponential scaling.

## 4 Using Modularity to Generalize in High Dimensions

So far, we have shown that a theoretical model treating NN learning as linear regression can closely model the generalization trends of actual NNs. In our model, the number of training samples required to generalize on a task with $m$ dimensional input scales exponentially with $m$. Now, we demonstrate that modular NNs (in contrast to the monolithic NNs studied so far) can avoid this exponential dependence on $m$ for tasks with an underlying modular structure. We will consider modular networks in which model parameters are divided into separate modules, each of which processes a projection of the input; monolithic (or nonmodular) networks in this context will correspond to networks without an explicit separation of parameters into modules. We first demonstrate the theoretical advantages of modular NNs under a specific form of modularity, then develop a modular NN learning rule to learn the underlying modular structure of a task. We then empirically validate our approach and demonstrate that our approach can learn the true modules underlying the task.

### 4.1 Sample Complexity of Modular Learning

Recall in Sec 3, modeling a NN as a linear function of its parameters successfully captured its generalization properties. We aim to use this model to demonstrate the improved generalization of modular learning. For analytical traceability, we restrict our analysis to a specific modular setting that captures crucial aspects of many real-world modular learning scenarios. In practical settings, modules often handle low-dimensional inputs, such as attention maps in Andreas et al. (). As such, we assume modules receive projected versions of the task input $x\in\mathbb{R}^{m}$. Our theoretical analysis will assume linear projections, but our method and experiments are also applied to non-linear projections. Moreover, we assume that the module outputs are summed to produce a final output, a feature of architectures such as Mixture of Experts (Shazeer et al., ). Specifically,
consider a modular NN constructed as a linear combination of general NNs (each constituting a module) of low-dimensional projections of the input:

$$ $\hat{y}(x)=\frac{1}{\sqrt{K}}\sum_{j=1}^{K}\hat{y}_{j}(\hat{U}_{j}^{T}x),$ (3) $$

where $\hat{U}_{j}\in\mathbb{R}^{m\times b}$ is a linear projection, $\hat{y}_{j}:\mathbb{R}^{b}\to\mathbb{R}^{d}$ is a NN. We will assume that $b$ is small, creating a bottleneck to each module’s input. We normalize by $\frac{1}{\sqrt{K}}$ to make the scale of $\hat{y}(x)$ invariant to $K$ (treating each term $\hat{y}_{j}(\hat{U}_{j}^{T}x)$ as independent, the sum of the terms has variance $O(K)$, so dividing my $\sqrt{K}$ makes their variance constant). Note that each module $\hat{y}_{j}$ is itself a monolithic NN with arbitary architecture; we do not restrict the form of the modules themselves. Using the linearity assumption of Eqn 18 (namely, the model is a linear function of parameters), we may model this as:

$$ $\hat{y}_{j}(\hat{U}_{j}^{T}x)=\varphi^{(U)}(x)\mathbf{F}(\hat{U}_{j})+\varphi^ {(W)}(x)\begin{bmatrix}I\\ 0\end{bmatrix}\theta_{j},$ (4) $$

where $\theta_{j}\in\mathbb{R}^{1\times p}$ are the parameters of $\hat{y}_{j}$, $\varphi^{(U)}(x)\in\mathbb{R}^{d\times mb}$ and $\varphi^{(W)}(x)\in\mathbb{R}^{d\times P}$ are feature matrices, $\mathbf{F}(\cdot)$ denotes flattening a matrix into a vector, and we consider the limit when $P\to\infty$. Observe that this is derived simply by assuming the model output is linear in $U_{j}$ and $\theta_{j}$ and defining the coefficients multiplying them as features $\phi^{(U)}(x)$ and $\phi^{(W)}(x)$. As before, we assume the features are distributed i.i.d from a unit Gaussian: $\varphi^{(U)}(x)_{i}\sim\mathcal{N}(0,I),\varphi^{(W)}(x)_{i}\sim\mathcal{N}(0
,I)$.
We may then write $\hat{y}$ as a linear model:

$$ $\hat{y}(x)=\varphi^{(U)}(x)\frac{1}{\sqrt{K}}\sum_{j=1}^{K}\mathbf{F}(\hat{U}_ {j})+\varphi^{(W)}(x)\begin{bmatrix}I\\ 0\end{bmatrix}\frac{1}{\sqrt{K}}\sum_{j=1}^{K}\theta_{j}.$ (5) $$

Next, we assume that the regression target $y$ has the same modular structure with $k$ modules:

$$ $y(x)=\frac{1}{\sqrt{k}}\sum_{j=1}^{k}y_{j}(U_{j}^{T}x),$ (6) $$

where $U_{j}\in\mathbb{R}^{m\times b}$ are the true projection directions and $y_{j}:\mathbb{R}^{b}\to\mathbb{R}^{c}$ are the true modules underlying the target. We apply the same linearity assumption to find:

$$ $y(x)=\varphi^{(U)}(x)\frac{1}{\sqrt{k}}\sum_{j=1}^{k}U_{j}+\varphi^{(W)}(x) \frac{1}{\sqrt{k}}\sum_{j=1}^{k}W_{j}$ (7) $$

where $W_{j}\in\mathbb{R}^{P}$ are the parameters of $y_{j}$. As in the case of monolithic networks, we assume $\mathbb{E}[W_{j}]=0$ and $\mathbb{E}[W_{j}W_{j}^{T}]=\Lambda$ where $\Lambda$ is diagonal. In this case, observe that each $W_{j}$ parameterizes a function with a $b$-dimensional input; thus it is appropriate to assume that $W_{j}$ is distributed as if $m=b$:

$$ $\lambda_{i}=c\left[i^{-\Omega^{-b}}-(i+1)^{-\Omega^{-b}}\right]$ (8) $$

To preserve spherical symmetry in the distribution of $U_{j}$, we assume that $\mathbb{E}[U_{j}]=0$ and $\mathbb{E}[\mathbf{F}(u_{j})\mathbf{F}(u_{j})^{T}]=I$. To complete our model definition, we define $\varphi(x)\in\mathbb{R}^{d\times(mb+P)}$ as the concatenation of $\varphi^{(U)}(x)$ and $\varphi^{(W)}(x)$:
$\varphi(x)=\left[\varphi^{(U)}(x),\varphi^{(W)}(x)\right]$.
We may then write:

$$ $y(x)=\varphi(x)\begin{bmatrix}I\\ 0\end{bmatrix}\overline{\theta}$ (9) $$

where $\overline{\theta}\in\mathbb{R}^{mb+p}$ is defined as: $\overline{\theta}=\frac{1}{\sqrt{K}}\left[\sum_{j=1}^{K}\mathbf{F}(\hat{U}_{j}
),\sum_{j=1}^{K}\theta_{j}\right]$. Similarly, we may write: $\hat{y}(x)=\varphi(x)\overline{W}$,
where $\overline{W}\in\mathbb{R}^{mb+P}$ is defined as:
$\overline{W}=\frac{1}{\sqrt{k}}\left[\sum_{j=1}^{k}\mathbf{F}(U_{j}),\sum_{j=1
}^{k}W_{j}\right]$.
Note that this model is nearly identical to that of monolithic NNs in Sec 3: the key difference is the different distribution of $\overline{W}$. $\overline{W}$ has covariance $\overline{\Lambda}=\begin{bmatrix}I&0\\
0&\Lambda\end{bmatrix}$ where $\Lambda$ is not dependent on $m$; $\Lambda$ is parameterized analogously to Section 3. Assuming $\hat{y}$ is trained to minimize squared loss on a training set, we may adapt Theorem 1 to this setting to compute the expected training and test loss of modular networks:

###### Theorem 2 .

Given a target function $y$ and model $\hat{y}$ estimated as described above, in the limit that $P\to\infty$, the expected test loss when averaging over $x$ and $\overline{W}$ is:

$$ $\lim_{P\to\infty}\mathbb{E}\left[||y(x)-\hat{y}(x)||^{2}\right]=d\operatorname {Tr}(\overline{\Lambda}_{2})F(dn,p)-d\frac{\min(dn,p)}{p}\operatorname{Tr}( \overline{\Lambda}_{1})+d\operatorname{Tr}(\overline{\Lambda})\\ =dF(dn,p)c(p+1)^{-\Omega^{-b}}-d\frac{\min(dn,p)}{p}\left(mb+c-c(p+1)^{-\Omega ^{-b}}\right)+dmb+dc$ (10) $$

The expected training loss is:

$$ $\lim_{P\to\infty}\frac{1}{n}\mathbb{E}\left[\left\|y(X)-\hat{y}(X)\right\|_{2} ^{2}\right]=\frac{dn-\min(dn,p)}{n}\operatorname{Tr}(\bar{\Lambda}_{2})\\ =\frac{dn-\min(dn,p)}{n}c(p+1)^{-\Omega^{-b}}$ (11) $$

with $F(n,p)$ defined as:

$$ $F(n,p)=\mathbb{E}\left[\left\|R^{\dagger}\right\|_{F}^{2}\right]$ (12) $$

where $R\in\mathbb{R}^{n\times p}$ has elements drawn i.i.d. from $\mathcal{N}(0,1)$.

The proof simply applies Theorem 1 with a different covariance matrix for $\hat{W}$; see App C for the full proof. When $Tr(\bar{\Lambda}_{2})$ is independent of $m$, we see that unlike the monolithic network, the training loss does not depend on $m$, and the dependence of the test loss on $m$ is linear. This is because the module inputs have effective dimensionality $b$ instead of $m$ due to the bottleneck caused by the module projections $\hat{U}_{j}$. Furthermore, in the underparameterized regime ($dn>p$), the test loss becomes $(dF(dn,p)+d)\operatorname{Tr}(\overline{\Lambda}_{2})$ which does not depend on $m$, implying that the $n$ required to reach a specific loss can be bounded by a function of only $p$ (assuming that some value of $p$ can achieve the desired loss). Thus, the sample complexity of modular NNs is independent of the task dimensionality. This dimension independence holds regardless of the parameterization of $\lambda_{i}$; the only condition is that the modules must have a dimension-independent input bottleneck (i.e. the covariance $\Lambda$ of module parameters must be independent of $m$). This result suggests that, unlike monolithic NNs, modular NNs can scale to high-dimensional, modular problems without requiring intractable amounts of data.

###### Theorem 2 .

### 4.2 Modular Learning Rule

Inspired by the improved theoretical generalizability of modular NNs, and the finding that modular architectures trained with gradient descent on a task often cannot exploit these efficiencies (Csordás et al., ; Mittal et al., ), we develop a modular learning rule that practically exhibits this advantage. Importantly, we now relax the assumption that module input projections are linear.

We consider modular regression tasks with targets constructed as follows:

$$ $y(x)=\sum_{j=1}^{k}y_{j}(x;U_{j})$ (13) $$

where $y_{j}$ are functions that depend on a potentially nonlinear projection of the input $x$ as represented by $y_{j}$ depending on both module projections $U_{j}$ and inputs $x$. Observe that this generalizes the linear projections considered before (in Eqn 6). Suppose we aim to model the target function by approximating the $U_{j}$ with $\hat{U}_{j}$ and the $y_{j}$ with $\hat{y}_{j}$ (parameterized as a neural network). We propose a kernel-based rule to learn the initializations of $\hat{U}_{j}$ from the training data; this allows us to efficiently learn the modules $\hat{y}_{j}$. Assume we are provided a set of training data $y(X)\in\mathbb{R}^{dn\times 1}$. Given the modular structure, we first aim to approximate the data as:

$$ $y(X)\approx\sum_{i=1}^{K}\varphi(X;\hat{U}_{i})\theta_{i},$ (14) $$

where $X\in\mathbb{R}^{n\times m}$, and $\varphi$ is an arbitrary nonlinearity applied elementwise to the input data such that $\varphi(X;\hat{U}_{i})\in\mathbb{R}^{dn\times p}$, $\theta_{i}\in\mathbb{R}^{p}$ and $K$ is the number of expected modules.

Figure: (a) Sine Wave Regression Task
Refer to caption: extracted/5841346/figs/compare_method1_sample_complexity1.5_edited.png

We expect that if $\hat{U}_{i}=U_{i}$ and $\varphi$ is sufficiently expressive, then $y(X)$ can be well approximated. Assuming $pK>dn$, observe that the minimum norm solution for $\theta_{i}$ can be computed as:

$$ $\begin{bmatrix}\theta_{1}\\ \theta_{2}\\ \vdots\\ \theta_{K}\end{bmatrix}=\begin{bmatrix}\varphi(X;\hat{U}_{1})&\varphi(X;\hat{U }_{2})&\cdots&\varphi(X;\hat{U}_{K})\end{bmatrix}^{\dagger}y(X)$ (15) $$

In general, such a solution exists for any choice of $\hat{U}_{i}$. However, we hypothesize that if the $\hat{U}_{i}$ is far from $U_{i}$, then the norm of the $\theta_{i}$ will be large: intuitively, interpolating the data along the ”incorrect” projection directions $\hat{U}_{i}$ will be more difficult. Thus, we optimize $\hat{U}_{i}$ to minimize the squared norm of $\theta_{i}$. Specifically, we minimize:

$$ $\sum_{i=1}^{K}\|\theta_{i}\|_{2}^{2}=y(X)^{T}\mathbf{K}^{-1}y(X)$ (16) $$

where $\mathbf{K}\in\mathbb{R}^{dn\times dn}$ is a kernel matrix applied to the data corresponding to the following kernel: $\mathbf{K}(x_{1},x_{2})=\sum_{i=1}^{K}\varphi(x_{1}^{T};\hat{U}_{i})\varphi(x_
{2}^{T};\hat{U}_{i})^{T}=\kappa(x_{1},x_{2};\hat{U}_{i})$,
where $\kappa(x_{1},x_{2};\hat{U}_{i})$ is a module-conditional kernel between $x_{1}$ and $x_{2}$. Experimentally, we tailor $\kappa$ to the modular structure of the problem we consider. Note that the above analysis only applies when $\phi$ is sufficiently expressive (i.e. $pK>dn$), which is a natural assumption for typically overparameterized models like neural networks. When models do not satisfy this assumption, Eqn 15 yields a solution minimizing the (generally nonzero) error between the predicted and true $y(X)$. Importantly, in this case, minimizing the norm of $\theta$ with respect to $\hat{U}_{i}$ may not necessarily yield a lower prediction error.

App E describes the specific choice of $\kappa$. Alg 1 shows the full procedure to find a single module projection $\hat{U}_{i}$; each step of the algorithm simply applies gradient on Eqn 16 with respect to $\hat{U}_{i}$. We repeat this procedure $K$ times with different random initializations to find the initial values of all $K$ module projections in our architecture. Then, we train all module parameters (including the $\hat{U}_{i}$) via gradient descent on the task loss. We stress that our approach is applicable to a fairly general set of modular architectures of the form $\sum_{j}\hat{y}_{j}(x;\hat{U}_{j})$: it does not restrict modules to receive only linear projections of inputs, and, moreover, does not restrict the form of the modules.

### 4.3 Experimental Results

Figure: Figure 5: Average cosine similarity between learned and target module directions over training for a modular NN initialized with our method vs. random initialization (baseline).
Refer to caption: extracted/5841346/figs/sims.png

We evaluate the generalizability of our method on a modular NN vs. baselines of a randomly initialized monolithic and modular NN trained on 1) sine wave regression tasks of varying dimensionality (fixing $k=m$), 2) a nonlinear variant of the sine wave regression task where the task has a nonlinear module structure, and 3) Compositional CIFAR-10 (based on Compositional MNIST (Jarvis et al., )), a modular task in which each input consists of multiple CIFAR-10 images and the goal is to simultaneously predict the classes of all images; see App E Fig 7 for an illustration. In Compositional CIFAR-10, each input is constructed as a concatenation of $k$ flattened CIFAR-10 images (resulting in a $3072k$ dimensional vector) and target outputs are $k$-hot encoded $10k$ dimensional vectors encoding the class of each component image. The modular architectures are constructed as $\hat{y}(x)=\sum_{j=1}^{k}\hat{y}_{j}(x;U_{j})$ where $y_{j}$ are fully connected, ReLU-activated networks, and the projection operation $U_{j}$ parameterizes the first layer. Monolithic architectures are normal fully connected ReLU-activated networks. See App E for further details on the datasets and the full experimental setup.

#### Modular NNs empirically generalize better in and out-of-distribution

As shown in Fig 4, our modular method generalizes better compared with both the monolithic baseline method and the modular baseline method as evaluated by sample complexity for the sine wave regression task and accuracy for Compositional CIFAR-10. On both tasks, our method’s advantage persists even on higher-dimensional inputs. Interestingly, on the sine-wave task, the monolithic baseline outperforms the modular baseline, highlighting the difficulty of optimizing modular architectures. We also conduct experiments on additional variants of Compositional CIFAR-10 that test out-of-distribution generalization: 1) our method learns to classify unseen class combinations, thus generalizing combinatorially (App F Tab 3) and 2) our method is robust to small amounts of Gaussian noise added to training inputs (App F Tab 4), thus generalizing to small distribution shifts.

#### Our learning rule finds the true task modules

Fig 5 computes a similarity score between our learned module projections ($\hat{U}$) and the target module projections ($U$) on the sine wave regression task; our method indeed aligns with the target modules. App F Fig 11 plots a low-dimensional representation of the target ($U$) and learned ($\hat{U}$) module projections: our learned NN initializations closely cluster around the target modules without any task training. On Compositional CIFAR-10, the learned module projections can be directly visualized as in App F Fig 12; here, our module initializations make each module sensitive to single component images without any task training, suggesting that our approach promotes generalization by correctly learning the modular structure of a task. We also conduct ablation studies in App F Fig 9 which show that 1) our method performs nearly as well as using the ground-truth module directions and 2) allowing our learned module directions to adapt directly to the task loss improves performance.

**Table 2: Comparison of our method with baselines on a nonlinear variant of the sine wave regression task where $k=m=5$. Test loss values are evaluated with standard errors over $5$ trials.**
| Method | Test Loss |
| --- | --- |
| Baseline monolithic | $1.302\pm 0.223$ |
| Baseline modular | $0.555\pm 0.136$ |
| Our method | $0.393\pm 0.099$ |

#### Our learning rule extends to nonlinear module projections

So far, we have considered modules ($y_{j}$ in Eqn 13) for which the input is a linear function of both $U_{j}$ and $x$. Next, we consider a nonlinear variant of the sine wave regression task in which modules are a function of $||u_{j}-x||_{2}$, where $u_{j}$ are the vector module projection directions, which is nonlinear in both $u_{j}$ and $x$. The modular architecture is constructed non-linearly as: $\frac{1}{\sqrt{K}}\sum_{j=1}^{K}\hat{y}_{j}(||\hat{u}_{j}-x||_{2})$; see App E for further details. Tab 2 illustrates that our method significantly outperforms the baselines, indicating that our method extends to nonlinear settings as well.

## 5 Discussion

Existing NN scaling laws show that in order to generalize on a task, monolithic NNs require an exponential number of training samples with the task’s dimensionality. In this paper, we develop theory demonstrating that modular NN can break this scaling law: they only require only a constant number of samples to generalize in terms of task dimension when applied to modular tasks. To our knowledge, we are the first to demonstrate such a result using explicit expressions for generalization error in modular NNs. Based on this theoretical finding, we propose a novel learning rule for modular NNs and demonstrate its improved generalization, both in and out of distribution, on a sine wave regression task and Compositional CIFAR-10.

In pursuit of explicit, non-asymptotic expressions for generalization error in modular and monolithic NNs, we make strong theoretical assumptions consistent with previous literature, which we hope can be further relaxed in future work. Moreover, our results apply to a specific form of modularity that captures structures in common real-world modular architectures but is not fully general. Notably, our theory considers linear module projections (although our method is also applied to nonlinear projections), and both the theory and the experiments assume that the model output is the sum of module outputs. We expect that future analyses can demonstrate the benefits of modularity more widely. For example, routing mechanisms are a popular type of modularity in which modules are flexibly composed or combined based on a routing network. We expect analyses similar to ours to show that, to the extent that routing mechanisms allow process low-dimensional inputs rather than the full task input, they also generalize better. Finally, we find that while the theory predicts a task-dimension-independent sample complexity for modular NNs, empirically we do not eliminate this dependence due to the difficulty of optimizing modular NNs in high dimensions. Nevertheless, our learning rule significantly eases this challenge.

Practically, we expect that modularity provides the most benefit for modular tasks with high-dimensional inputs; this is because the relative sample complexity improvement between nonmodular and modular tasks is greater when task dimensionality increases. Indeed, as discussed in Sec 2, modularity empirically significantly improves generalization in domains ranging from reinforcement learning and robotics to visual question answering and language modeling, all of which can be highly compositional and can have high-dimensional task inputs. We speculate that the modular structure of self-attention-based architectures may explain their success in many of these domains. Our findings provide a step toward fundamentally understanding how modularity can be better applied to solve high-dimensional generalization problems.

## Appendix A Theoretical Model of Neural Network Generalization

### A.1 Setup

We consider a regression task with input $x\in\mathbb{R}^{m}$ and a feature matrix $\varphi(x)\in\mathbb{R}^{d\times P}$ such that over the data distribution, the features are distributed i.i.d. from a unit Gaussian: $\varphi(x)_{i,j}\sim\mathcal{N}(0,1)$. We consider the limit when $P\to\infty$. Suppose that our regression target function $y:\mathbb{R}^{m}\to\mathbb{R}^{d}$ is constructed linearly from $\varphi(x)$:

$$ $y(x)=\varphi(x)W,$ (17) $$

where $W\in\mathbb{R}^{P\times 1}$. To accommodate multidimensional outputs, note that we shape input features $\phi(x)$ as a matrix (with one row for each output dimension) and parameters $W$ as a vector. This choice is more general than parameterizing each output dimension independently (this can be captured as a special case of our approach) and , moreover, aligns with prior theoretical literature (Jacot et al., ). Assume $\mathbb{E}[W]=0$ and $\mathbb{E}[WW^{T}]=\Lambda$, where $\Lambda$ is diagonal and $\operatorname{Tr}(\Lambda)$ is finite. Suppose we aim to approximate $y(x)$ using a model $\hat{y}(x)$ constructed as follows:

$$ $\hat{y}(x)=\varphi(x)\begin{bmatrix}I\\ 0\end{bmatrix}\theta,$ (18) $$

where $\theta\in\mathbb{R}^{p\times 1}$ are model parameters and $p$ is the number of parameters. This corresponds to the model only being able to control $p$ of the $P$ true underlying parameters in the construction of $y$. We decompose $\Lambda$ blockwise as: $\Lambda=\begin{bmatrix}\Lambda_{1}&0\\
0&\Lambda_{2}\end{bmatrix}$
where $\Lambda_{1}\in\mathbb{R}^{p\times p}$ and $\Lambda_{2}\in\mathbb{R}^{(P-p)\times(P-p)}$.

We make a specific choice of parameterization for the individual elements $\lambda_{i}$ of $\Lambda$:

$$ $\lambda_{i}=c\left[i^{-\Omega^{-m}}-(i+1)^{-\Omega^{-m}}\right]$ (19) $$

for some constants $c$ and $\Omega$. We justify this choice as follows: we define the effective dimensionality of $\Lambda$ as $\frac{(\sum_{i}\lambda_{i})^{2}}{\sum_{i}\lambda_{i}^{2}}$ (this measure approximates the $\ell_{0}$ norm (Krishnan et al., ), and thus can be used as a measure of $\Lambda$’s dimensionality). For large $m$, this can be approximated as: $\frac{(\sum_{i}\lambda_{i})^{2}}{\sum_{i}\lambda_{i}^{2}}\approx\Omega^{2m}$; we interpret this as $y$ having $\Omega^{2m}$ free parameters. This is consistent with the observation that regression on an $m$-dimensional input space has a function space that scales exponentially with $m$ (McRae et al., ). Intuitively, this is because the function must have enough free parameters to express values at all points in volume of its input space, and volume scales exponentially with $m$ (Sharma & Kaplan, ).

Next, suppose we are given a set of training data $X$ with associated feature matrix transformation $\varphi(X)\in\mathbb{R}^{dn\times P}$, where $n$ is the number of data points. Suppose $\theta$ is optimized to find the minimum norm interpolating solution to the data:

$$ $\min_{\theta}\left\{\left\|\varphi(X)\begin{bmatrix}I\\ 0\end{bmatrix}\theta-\varphi(X)W\right\|_{2}^{2}+\frac{\gamma}{2}\|\theta\|_{2 }^{2}\right\},$ (20) $$

where $\gamma\to 0$. Recall that the solution can be found as:

$$ $\theta=\left[\varphi(X)\begin{bmatrix}I\\ 0\end{bmatrix}\right]^{\dagger}\varphi(X)W.$ (21) $$

Thus, the model’s prediction on a point $x$ is:

$$ $\hat{y}(x)=\varphi(x)\begin{bmatrix}I\\ 0\end{bmatrix}\left[\varphi(X)\begin{bmatrix}I\\ 0\end{bmatrix}\right]^{\dagger}\varphi(X)W.$ (22) $$

### A.2 Empirical Validation

Next, we empirically validate our theoretical model on a parametrically-variable modular sine wave regression task (see App E for task details). The task allows us to quantify how NN generalization depends on the number of dimensions of variation such as the dimensionality $m$ of the task input, the number of model parameters $p$, the number of samples $n$ and the number of modules $k$ in the construction of the target function. Note that our theoretical model does not include the number of modules $k$ since it does not explicitly construct the target modularly. Thus, directly applying our theory, we would expect the loss to be invariant to $k$. Intuitively, this is because varying $k$ increases the complexity of the task in a way that is irrelevant for generalization.

#### Trends of generalization error of NNs

Now, we train NNs of various architectures on the task and observe error trends as a function of $k,m,p$ and $n$. We fit parameters $c$ and $\Omega$ of our theoretical model to our task. Furthermore, because each parameter in our NNs may not correspond to a single parameter in our theoretical model, we use a linear scaling of the number of true NN parameters to estimate the number of effective parameters in our theoretical model; specifically, we estimate $p=\alpha p^{\prime}$ where $p^{\prime}$ is the actual number of NN parameters. See App E for more details.

Fig 2 shows that our theoretical model can capture many trends of the training and test loss as a function of $k,m,p$ and $n$. In particular, our model predicts the invariance of loss to $k$, the sub-linear increase in loss with $m$, and the double descent behavior of loss with $p$ and $n$. Notably, our model predicts the empirical location of the interpolation threshold as seen in the last two plots of Fig 2.

We note two key discrepancies between our theory and empirical results: first, the loss is empirically larger than predicted for small amounts of training data, and second, the error spike at the interpolation threshold is smaller than predicted by the theory.

We believe the first discrepancy is due to imperfect optimization of neural networks, especially in low data regimes. Note that the linearized analysis assumes that the linear model solution finds the exact global optimum. However, the actual optimization landscape for modular architectures is highly non-convex, and the global optimum may not be found especially for small datasets (indeed, we find a significant discrepancy between predicted and actual training loss values for small data size $n$; in the overparameterized regime, the predicted training error is exactly 0 0). We believe this causes the discrepancy between predicted and actual test errors in low data regimes.

We hypothesize that the second discrepancy is also partly due to imperfect optimization. This is because the interpolation threshold spike can be viewed as highly adverse fitting to spurious training set patterns. This imperfect optimization is more pronounced at smaller $m$. Despite these discrepancies, we nevertheless find that our theory precisely captures the key trends of empirical test error.

Finally, we consider the trend between $m$ and $n$ implied by our model. As Fig 3 reveals, for various nonmodular architectures, the sample complexity grows approximately exponentially with the task dimensionality, consistent with prior theoretical observations (see Sec 2.2). This implies that generalizing on high-dimensional problems can require a massive number of samples.

## Appendix B Proof of Theorem 1

###### Proof.

Test set error We first compute the expected test set error. Note that the squared error can be written as:

$$ $||\hat{y}(x)-y(x)||^{2}=\left|\left|\varphi(x)\begin{bmatrix}I\\ 0\end{bmatrix}\left[\varphi(X)\begin{bmatrix}I\\ 0\end{bmatrix}\right]^{\dagger}\varphi(X)W-\varphi(x)W\right|\right|^{2}\\ =\left|\left|\varphi(x)\left[\begin{bmatrix}I\\ 0\end{bmatrix}\left[\varphi(X)\begin{bmatrix}I\\ 0\end{bmatrix}\right]^{\dagger}\varphi(X)-I\right]W\right|\right|^{2}$ (23) $$

For notational convenience, define $A$ as the first $p$ columns of $\varphi(X)$ and $B$ as the remaining $P-p$ columns such that $\varphi(X)=[A,B]$. Also, define $M=\begin{bmatrix}I\\
0\end{bmatrix}\left[\varphi(X)\begin{bmatrix}I\\
0\end{bmatrix}\right]^{\dagger}\varphi(X)=\begin{bmatrix}I\\
0\end{bmatrix}A^{\dagger}[A,B]=\begin{bmatrix}A^{\dagger}A&A^{\dagger}B\\
0&0\end{bmatrix}$. Then, using the cyclic property of trace, the squared error can be written as

$$ $||\varphi(x)(M-I)W||^{2}=\operatorname{Tr}\left(\varphi(x)(M-I)WW^{T}(M-I)^{T} \varphi(x)^{T}\right)=\operatorname{Tr}\left((M-I)WW^{T}(M-I)^{T}\varphi(x)^{T }\varphi(x)\right)$ (24) $$

Next, we can take the expectation with to $x$ and $W$ and use the fact that $\mathbb{E}[\varphi(x)^{T}\varphi(x)]=dI$ and $\mathbb{E}[WW^{T}]=\Lambda$ to find that:

$$ $\mathbb{E}\left[||\hat{y}(x)-y(x)||^{2}\right]=d\mathbb{E}\left[\operatorname{ Tr}((M-I)\Lambda(M-I)^{T})\right]$ (25) $$

Finally, expanding:

$$ $d\mathbb{E}\left[\operatorname{Tr}((M-I)\Lambda(M-I)^{T})\right]=d\mathbb{E} \left[\operatorname{Tr}(M^{T}M\Lambda)\right]-2d\mathbb{E}\left[\operatorname{ Tr}(M\Lambda)\right]+d\operatorname{Tr}(\Lambda)$ (26) $$

Next, we compute $\mathbb{E}\left[\operatorname{Tr}(M^{T}M\Lambda)\right]$:

$$ $\mathbb{E}\left[\operatorname{Tr}(M^{T}M\Lambda)\right]=\mathbb{E}\left[ \operatorname{Tr}\left(\begin{bmatrix}A^{\dagger}A&0\\ B^{T}A^{\dagger T}&0\end{bmatrix}\begin{bmatrix}A^{\dagger}A&A^{\dagger}B\\ 0&0\end{bmatrix}\Lambda\right)\right]\\ =\mathbb{E}\left[\operatorname{Tr}\left(\begin{bmatrix}A^{\dagger}A&A^{\dagger }AA^{\dagger}B\\ B^{T}A^{\dagger T}A^{\dagger}A&B^{T}A^{\dagger T}A^{\dagger}B\end{bmatrix} \Lambda\right)\right]$ (27) $$

Next, we decompose $\Lambda$ blockwise as:

$$ $\Lambda=\begin{bmatrix}\Lambda_{1}&0\\ 0&\Lambda_{2}\end{bmatrix}$ (28) $$

where $\Lambda_{1}\in\mathbb{R}^{p\times p}$ and $\Lambda_{2}\in\mathbb{R}^{(P-p)\times(P-p)}$. Then:

$$ $\mathbb{E}\left[\operatorname{Tr}(M^{T}M\Lambda)\right]=\mathbb{E}\left[ \operatorname{Tr}\left(\begin{bmatrix}A^{\dagger}A\Lambda_{1}&A^{\dagger}AA^{ \dagger}B\Lambda_{2}\\ B^{T}A^{\dagger T}A^{\dagger}A\Lambda_{1}&B^{T}A^{\dagger T}A^{\dagger}B \Lambda_{2}\end{bmatrix}\right)\right]\\ =\mathbb{E}\left[\operatorname{Tr}(A^{\dagger}A\Lambda_{1})\right]+\mathbb{E} \left[\operatorname{Tr}(A^{\dagger T}A^{\dagger}B\Lambda_{2}B^{T})\right]$ (29) $$

Next, consider, $\mathbb{E}[\operatorname{Tr}(M\Lambda)]$:

$$ $\mathbb{E}[\operatorname{Tr}(M\Lambda)]=\mathbb{E}\left[\operatorname{Tr}\left (\begin{bmatrix}A^{\dagger}A&A^{\dagger}B\\ 0&0\end{bmatrix}\Lambda\right)\right]=\mathbb{E}\left[\operatorname{Tr}(A^{ \dagger}A\Lambda_{1})\right]$ (30) $$

Combining this result with the earlier result, the expected squared error can be expressed as:

$$ $\mathbb{E}\left[||\hat{y}(x)-y(x)||^{2}\right]=d\mathbb{E}\left[\operatorname{ Tr}(A^{\dagger T}A^{\dagger}B\Lambda_{2}B^{T})\right]-d\mathbb{E}\left[ \operatorname{Tr}(A^{\dagger}A\Lambda_{1})\right]+d\operatorname{Tr}(\Lambda)$ (31) $$

Next, we evaluate $\mathbb{E}\left[\operatorname{Tr}(A^{\dagger T}A^{\dagger}B\Lambda_{2}B^{T})\right]$. By linearity of trace, and the independence of $A$ and $B$, we have:

$$ $\mathbb{E}\left[\operatorname{Tr}(A^{\dagger T}A^{\dagger}B\Lambda_{2}B^{T}) \right]=\operatorname{Tr}\left(\mathbb{E}[A^{\dagger T}A^{\dagger}]\cdot \mathbb{E}[B\Lambda_{2}B^{T}]\right)$ (32) $$

Define the following quantities:

$$ $\alpha=\mathbb{E}\left[A^{\dagger T}_{:,1}A^{\dagger}_{:,1}\right]$ (33) $$

and

$$ $\beta=\mathbb{E}\left[B_{1,:}\Lambda_{2}B^{T}_{1,:}\right]$ (34) $$

where $\cdot_{i,j}$ indicates the $i$th row and $j$th column of the argument. Note that by symmetry over the data points and output dimensions, both $\mathbb{E}\left[A^{\dagger T}A^{\dagger}\right]\in\mathbb{R}^{dn\times dn}$ and $\mathbb{E}\left[B\Lambda_{2}B^{T}\right]\in\mathbb{R}^{dn\times dn}$ must be proportional to the identity matrix. Thus, the expectation of their top left entry is the same as the expectation of any other entry:

$$ $\mathbb{E}\left[A^{\dagger T}A^{\dagger}\right]=\alpha I$ (35) $$

and

$$ $\mathbb{E}\left[B\Lambda_{2}B^{T}\right]=\beta I$ (36) $$

Then,

$$ $\mathbb{E}\left[\operatorname{Tr}(A^{\dagger T}A^{\dagger}B\Lambda_{2}B^{T}) \right]=\operatorname{Tr}(\alpha I\beta I)=\alpha\beta n$ (37) $$

To evaluate $\alpha$, observe that since $A$ has elements distributed from $\mathcal{N}(0,1)$:

$$ $F(dn,p)=\mathbb{E}\left[\|A^{\dagger}\|_{F}^{2}\right]=\operatorname{Tr}\left( \mathbb{E}[A^{\dagger T}A^{\dagger}]\right)$ (38) $$

Using the definition of $\alpha$:

$$ $F(dn,p)=\alpha dn$ (39) $$

Therefore, $\alpha=\frac{F(dn,p)}{dn}$. To evaluate $\beta$, note that $B_{1,:}$ has elements distributed from $\mathcal{N}(0,1)$. Thus, $\beta$ is simply:

$$ $\beta=\mathbb{E}\left[B_{1,:}\Lambda_{2}B^{T}_{1,:}\right]=\sum_{i=1}^{P-p} \mathbb{E}[B_{1,i}^{2}]\Lambda_{2,i}=\operatorname{Tr}(\Lambda_{2})$ (40) $$

Substituting $\alpha$ and $\beta$ into the expression for $\mathbb{E}\left[\operatorname{Tr}(A^{\dagger T}A^{\dagger}B\Lambda_{2}B^{T})\right]$:

$$ $\mathbb{E}\left[\operatorname{Tr}(A^{\dagger T}A^{\dagger}B\Lambda_{2}B^{T}) \right]=\operatorname{Tr}(\Lambda_{2})F(dn,p)$ (41) $$

Next, we evaluate $\mathbb{E}\left[\operatorname{Tr}(A^{\dagger}A\Lambda_{1})\right]$. First, define the singular value decomposition of $A$ as $A=U\Sigma V^{T}$. Then, expanding $A$ and using the cyclic property of trace:

$$ $\mathbb{E}\left[\operatorname{Tr}(A^{\dagger}A\Lambda_{1})\right]=\mathbb{E} \left[\operatorname{Tr}(V\Sigma^{\dagger}U^{T}U\Sigma V^{T}\Lambda_{1})\right] =\mathbb{E}\left[\operatorname{Tr}(V\Sigma^{\dagger}\Sigma V^{T}\Lambda_{1})\right]$ (42) $$

Applying the linearity of trace:

$$ $\mathbb{E}\left[\operatorname{Tr}(A^{\dagger}A\Lambda_{1})\right]= \operatorname{Tr}\left(\mathbb{E}[V\Sigma^{\dagger}\Sigma V^{T}]\Lambda_{1}\right)$ (43) $$

Now, we examine $\mathbb{E}\left[V\Sigma^{\dagger}\Sigma V^{T}\right]\in\mathbb{R}^{p\times p}$. First, note that $\Sigma^{\dagger}\Sigma$ is a diagonal matrix with entries $1$ and 0 0: specifically, it has $\min(dn,p)$ $1$s and remaining entries (if any) 0 0. Thus, we may write:

$$ $\mathbb{E}\left[V\Sigma^{\dagger}\Sigma V^{T}\right]=\sum_{i=1}^{\min(dn,p)} \mathbb{E}\left[V_{:,i}V_{:,i}^{T}\right]$ (44) $$

Next, note that the distribution of $A$ is symmetric to rotations of its $p$ columns. Thus $V_{:,i}$ must also have a rotationally symmetric distribution. Since $\|V_{:,i}\|_{2}=1$, $\mathbb{E}\left[V_{:,i}V_{:,i}^{T}\right]=\frac{1}{p}I$. Thus:

$$ $\mathbb{E}\left[V\Sigma^{\dagger}\Sigma V^{T}\right]=\frac{\min(dn,p)}{p}I$ (45) $$

Substituting into the expression for $\mathbb{E}\left[\operatorname{Tr}(A^{\dagger}A\Lambda_{1})\right]$:

$$ $\mathbb{E}\left[\operatorname{Tr}(A^{\dagger}A\Lambda_{1})\right]=\frac{\min( dn,p)}{p}\operatorname{Tr}(\Lambda_{1})$ (46) $$

Combining the results from earlier, the expected squared error can be written as:

$$ $\mathbb{E}\left[||\hat{y}(x)-y(x)||^{2}\right]=d\operatorname{Tr}(\Lambda_{2}) F(dn,p)-d\frac{\min(dn,p)}{p}\operatorname{Tr}(\Lambda_{1})+d\operatorname{Tr} (\Lambda)$ (47) $$

Next, using the definition of $\lambda_{i}=c\left[i^{-\Omega^{-m}}-(i+1)^{-\Omega^{-m}}\right]$, observe that:

$$ $\operatorname{Tr}(\Lambda_{1})=\sum_{i=1}^{p}\lambda_{i}=\sum_{i=1}^{p}c\left[ i^{-\Omega^{-m}}-(i+1)^{-\Omega^{-m}}\right]=c-c(p+1)^{-\Omega^{-m}}$ (48) $$

$$ $\operatorname{Tr}(\Lambda_{2})=\sum_{i=p+1}^{\infty}\lambda_{i}=\sum_{i=p+1}^{ \infty}c\left[i^{-\Omega^{-m}}-(i+1)^{-\Omega^{-m}}\right]=c(p+1)^{-\Omega^{-m}}$ (49) $$

$$ $\operatorname{Tr}(\Lambda)=\operatorname{Tr}(\Lambda_{1})+\operatorname{Tr}( \Lambda_{2})=c$ (50) $$

We finally use the expressions for $\Lambda_{1}$ and $\Lambda_{2}$ to write the result in terms of $c$ and $\Omega$:

$$ $\mathbb{E}\left[||\hat{y}(x)-y(x)||^{2}\right]=dF(dn,p)c(p+1)^{-\Omega^{-m}}-d \frac{\min(dn,p)}{p}\left(c-c(p+1)^{-\Omega^{-m}}\right)+dc$ (51) $$

Training set error Now, we compute the training set error. Writing out the summed training set error over all data points:

$$ $\|\hat{y}(X)-y(X)\|_{2}^{2}=\left\|\varphi(X)\begin{bmatrix}I\\ 0\end{bmatrix}\left[\varphi(X)\begin{bmatrix}I\\ 0\end{bmatrix}\right]^{\dagger}\varphi(X)W-\varphi(X)W\right\|_{2}^{2}=\left\| (AA^{\dagger}-I)[A,B]W\right\|_{2}^{2}\\ =\left\|[(AA^{\dagger}-I)A,(AA^{\dagger}-I)B]W\right\|_{2}^{2}=\left\|[0,(AA^{ \dagger}-I)B]W\right\|_{2}^{2}$ (52) $$

Expressing this squared norm as a trace and using the cyclic property of trace:

$$ $\left\|[0,(AA^{\dagger}-I)B]W\right\|_{2}^{2}=\operatorname{Tr}\left(W^{T}[0,( AA^{\dagger}-I)B]^{T}[0,(AA^{\dagger}-I)B]W\right)\\ =\operatorname{Tr}\left([0,(AA^{\dagger}-I)B]^{T}[0,(AA^{\dagger}-I)B]WW^{T}\right)$ (53) $$

Taking the expectation with respect to $W$:

$$ $\mathbb{E}\left[\operatorname{Tr}\left([0,(AA^{\dagger}-I)B]^{T}[0,(AA^{ \dagger}-I)B]WW^{T}\right)\right]=\operatorname{Tr}\left(\mathbb{E}\left[B^{T} (AA^{\dagger}-I)^{2}B\right]\Lambda_{2}\right)$ (54) $$

Note that $(AA^{\dagger}-I)^{2}=I-AA^{\dagger}$
Again using the cyclic property of trace and the independence of $A$ and $B$, we find:

$$ $\operatorname{Tr}\left(\mathbb{E}[B^{T}(AA^{\dagger}-I)^{2}B]\Lambda_{2}\right )=\operatorname{Tr}\left(\mathbb{E}[I-AA^{\dagger}]\mathbb{E}[B\Lambda_{2}B^{T }]\right)$ (55) $$

From the calculations for test set error, we have:

$$ $\mathbb{E}\left[B\Lambda_{2}B^{T}\right]=\operatorname{Tr}(\Lambda_{2})I$ (56) $$

Substituting into the earlier expression:

$$ $\operatorname{Tr}\left(\mathbb{E}[I-AA^{\dagger}]\mathbb{E}[B\Lambda_{2}B^{T}] \right)=\operatorname{Tr}\left(\mathbb{E}[I-AA^{\dagger}]\operatorname{Tr}( \Lambda_{2})I\right)=\operatorname{Tr}(\Lambda_{2})\mathbb{E}[\operatorname{Tr }(I-AA^{\dagger})]$ (57) $$

To evaluate $\mathbb{E}[\operatorname{Tr}(I-AA^{\dagger})]$, we use the singular value decomposition of $A=U\Sigma V^{T}$ and the cyclic property of trace:

$$ $\mathbb{E}\left[\operatorname{Tr}(I-AA^{\dagger})\right]=\mathbb{E}\left[ \operatorname{Tr}(I-U\Sigma V^{T}V\Sigma^{\dagger}U^{T})\right]=\mathbb{E} \left[\operatorname{Tr}(I-U\Sigma\Sigma^{\dagger}U^{T})\right]=\mathbb{E}\left [\operatorname{Tr}(I-\Sigma\Sigma^{\dagger})\right]$ (58) $$

Observe that $A$ has full rank with probability $1$. Thus, $\Sigma$ also has full rank with probability $1$, implying that $\Sigma\Sigma^{\dagger}$ is with probability $1$ a diagonal matrix with $\min(dn,p)$ $1$s and remaining entries (if any) 0 0:

$$ $\mathbb{E}\left[\operatorname{Tr}(I-\Sigma\Sigma^{\dagger})\right]=dn-\min(dn,p)$ (59) $$

Substituting into the earlier expression, we have a result for the total training set error over all training points:

$$ $\mathbb{E}\left[\|\hat{y}(X)-y(X)\|_{2}^{2}\right]=\operatorname{Tr}(\Lambda_{ 2})(dn-\min(dn,p))$ (60) $$

To arrive at the final result for expected training set error we simply divide by $n$ and express $\operatorname{Tr}(\Lambda_{2})$ in terms of $c$ and $\Omega$

$$ $\frac{1}{n}\mathbb{E}[||\hat{y}(X)-y(X)||_{2}^{2}]=\frac{dn-\min(dn,p)}{n} \operatorname{Tr}(\Lambda_{2})=\frac{dn-\min(dn,p)}{n}c(p+1)^{-\Omega^{-m}}$ (61) $$

∎

###### Proof.

## Appendix C Proof of Theorem 2

###### Proof.

Test set error
Using the same techniques as in the proof of Theorem 1, we may write the expected test set error in terms of $\Lambda$ as:

$$ $\mathbb{E}\left[||\hat{y}(x)-y(x)||^{2}\right]=d\operatorname{Tr}(\overline{ \Lambda}_{2})F(dn,p)-d\frac{\min(dn,p)}{p}\operatorname{Tr}(\overline{\Lambda} _{1})+d\operatorname{Tr}(\overline{\Lambda})$ (62) $$

Using the definition $\lambda_{i}=c\left[i^{-\Omega^{-b}}-(i+1)^{-\Omega^{-b}}\right]$ and the assumption $p>mb$, observe that:

$$ $\operatorname{Tr}(\overline{\Lambda}_{1})=mb+\sum_{i=1}^{p}\lambda_{i}=mb+\sum _{i=1}^{p}c\left[i^{-\Omega^{-b}}-(i+1)^{-\Omega^{-b}}\right]=mb+c-c(p+1)^{- \Omega^{-b}}$ (63) $$

$$ $\operatorname{Tr}(\overline{\Lambda}_{2})=\sum_{i=p+1}^{\infty}\lambda_{i}= \sum_{i=p+1}^{\infty}c\left[i^{-\Omega^{-b}}-(i+1)^{-\Omega^{-b}}\right]=c(p+1 )^{-\Omega^{-b}}$ (64) $$

$$ $\operatorname{Tr}(\bar{\Lambda})=\operatorname{Tr}(\bar{\Lambda}_{1})+ \operatorname{Tr}(\bar{\Lambda}_{2})=mb+c$ (65) $$

Substituting these expressions, the expected test set error is:

$$ $\mathbb{E}\left[||\hat{y}(x)-y(x)||^{2}\right]=dF(dn,p)c(p+1)^{-\Omega^{-b}}-d \frac{\min(dn,p)}{p}\left(mb+c-c(p+1)^{-\Omega^{-b}}\right)+dmb+dc$ (66) $$

Training set error
Again, using the techniques in the proof of Theorem 1, we write the expected training set error in terms of $\bar{\Lambda}_{2}$

$$ $\frac{1}{n}\mathbb{E}\left[\|\hat{y}(X)-y(X)\|_{2}^{2}\right]=\frac{dn-\min(dn ,p)}{n}\operatorname{Tr}(\bar{\Lambda}_{2})$ (67) $$

Using the expression for $\operatorname{Tr}(\bar{\Lambda}_{2})$:

$$ $\frac{1}{n}\mathbb{E}\left[\|\hat{y}(X)-y(X)\|_{2}^{2}\right]=\frac{dn-\min(dn ,p)}{n}c(p+1)^{-\Omega^{-b}}$ (68) $$

∎

###### Proof.

## Appendix D Properties of F ⁢ ( n , p ) 𝐹 𝑛 𝑝 F(n,p) italic_F ( italic_n , italic_p )

In this section, we summarize some known properties about the function $F(n,p)$, which appears in Theorem 1. Recall that

$$ $F(n,p)=\mathbb{E}\left[\|R^{\dagger}\|_{F}^{2}\right],$ $$

where $R\in\mathbb{R}^{n\times p}$ has elements drawn i.i.d. from $\mathcal{N}(0,1)$.

In the regime $|n-p|\geq 2$, an exact closed form is given by

$$ $F(n,p)=\frac{\min(n,p)}{|n-p|-1}\qquad\text{if }|n-p|\geq 2.$ $$

Computing the square of the Frobenius norm of $R^{\dagger}$ is equivalent to finding $\operatorname{Tr}(R^{\dagger}R^{\dagger T})=\operatorname{Tr}(R^{\dagger T}R^{
\dagger})$.

When $n-p\geq 2$, $R^{\dagger}R^{\dagger T}$ is a $p\times p$ matrix with Inverse-Wishart distribution of identity covariance, which has mean $\frac{1}{n-p-1}I$ (Von Rosen, ). Therefore, the expected value of its trace is $\frac{n}{n-p-1}$. Analogously, when $p-n\geq 2$, $R^{\dagger T}R^{\dagger}$ is a $n\times n$ matrix with Inverse-Wishart distribution of identity covariance, which has mean $\frac{1}{p-n-1}I$, so the expected value of its trace is $\frac{n}{p-n-1}$.

In the case where $p=n$, bound on the Frobenius norm of $R^{\dagger}$ are known (Szarek, ). However, for the cases where $|p-n|\leq 1$, $F(n,p)$ has no known explicit form, so it was computed by averaging over 100 random trials.

## Appendix E Experimental Details

### E.1 Dataset Details

#### Sine Wave Regression Task

We construct a regression problem where the regression target is constructed as a sum of $k$ functions of one-dimensional linear projections of the input $x\in\mathbb{R}^{m}$. Specifically, the regression target function $y:\mathbb{R}^{m}\to\mathbb{R}$ is constructed as follows:

$$ $y(x)=\frac{1}{\sqrt{k}}\sum_{i=1}^{k}\sum_{j=1}^{\tau}a_{ij}\sin(\omega_{ij}u_ {i}^{T}x+\phi_{ij})$ (69) $$

where $a_{ij},\omega_{ij},\phi_{ij}\in\mathbb{R}$, $u_{i}\in\mathbb{R}^{m}$, $\tau$ is the number of sine waves that comprise each function of one-dimensional linear projection $u_{i}^{T}x$. We sample these parameters of the target function independently from the following distributions: $a_{ij}\sim\mathcal{N}(0,1)$, $\omega_{ij}\sim\mathcal{N}(0,4\pi^{2})$, $\phi_{ij}\sim\mathcal{U}(-\pi,\pi)$ where $\mathcal{N}$ denotes a Gaussian distribution and $\mathcal{U}$ denotes a uniform distribution. $u_{i}$ is drawn uniformly from a unit sphere. Note that the full function is made of $k$ separate functions of one-dimensional projections of $x$, each of which has $\tau$ sine components. We normalize by $\frac{1}{\sqrt{k}}$ so that $y(x)$ does not scale with $k$. Note that the task is modular: although the task takes an $m$ dimensional input, the target is constructed as a combination of $k$ modules operating on $1$ dimensional projections of the input. Given the projections, parameterizing functions of the projections are sufficient to parameterize the full task.

A training dataset is generated by first drawing $n$ training samples $x$ from a mean-zero Gaussian: $x\sim N(0,I)$. Then, for each $x$, the regression target $y(x)$ is computed. The test dataset is constructed analogously. See Fig 6 for an illustrated example target function in one dimension.

Recall that any square integrable function can be approximated on a finite interval to arbitrary precision by a sufficiently large Fourier series. Thus, we may expect that as $T$ approaches infinity, the functional form in Equation 69 can express any function constructed as a sum of square integrable functions of the $u_{i}^{T}x$: $\sum_{i=1}^{k}y_{i}(u_{i}^{T}x)$ for any square integrable $y_{i}$.

Note that $m$ controls the dimensionality of $x$; thus, we may make the task test generalization over arbitrarily high dimensions by simply increasing $m$. This is significant because prior work shows that the number of samples required to generalize to a fixed precision on a regression problem scales exponentially with the intrinsic dimensionality of the task input (McRae et al., ; Sharma & Kaplan, ). Thus, even with a relatively simple task construction, we may expect to produce tasks with arbitrary difficulty as measured by sample complexity.

Figure: Figure 6: An illustration of a single one-dimensional projection in our sine wave regression task with sampled training points. The target function in this projection direction is made of three sine wave components summed together.
Refer to caption: extracted/5841346/figs/sine_wave_diagram_edited.png

#### Nonlinear Sine Wave Regression Task

We also test our approach on a non-linear variation of our sine wave regression task. Recall that in the original sine wave regression task, outputs are constructed as:

$$ $y(x)=\frac{1}{\sqrt{k}}\sum_{i=1}^{k}\sum_{j=1}^{\tau}a_{ij}\sin(\omega_{ij}u_ {i}^{T}x+\phi_{ij})$ (70) $$

where $u_{i}$ are module projection directions. Note that this task has linear module input projections (the projections $u_{i}^{T}x$ are linear functions of $x$ and $u_{i}$). In our nonlinear variant, we consider the following outputs constructed with non-linear module input projections:

$$ $y(x)=\frac{1}{\sqrt{k}}\sum_{i=1}^{k}\sum_{j=1}^{\tau}a_{ij}\sin(\omega_{ij}|| u_{i}-x||_{2}+\phi_{ij})$ (71) $$

where $u_{i}^{T}x$ is replaced with $||u_{i}-x||_{2}$, which is non-linear in both $u_{i}$ and $x$. Remaining task parameters are set the same way as in the original sine wave regression task.

#### Compositional CIFAR-10

We conduct experiments on a Compositional CIFAR-10 dataset inspired by the Compositional MNIST dataset of (Jarvis et al., ). In the task, combinations of $k$ CIFAR-10 images are concatenated together and the model is asked to predict the class of all component images simultaneously. Fig 7 illustrates an example input. Inputs are flattened to remove all spatial structure. Outputs are $k$-hot encoded vectors constructed by concatenating the $1$-hot encoded labels for each component image; thus, targets are $10k$ dimensional. For this task, accuracies are reported on average over component images (for instance, if a model correctly guesses the class of two out of four images, the accuracy would be $50\%$).

Training and test sets for this task are constructed respectively as follows: each input in the training (or test) set is produced by randomly selecting $k$ images without replacement from the original CIFAR-10 training (or test) set and concatenating them in a random permutation. With $k$ images, there are $10^{k}$ possible class permutations for each input. We use a fixed training set size of $10^{6}$; thus, the probability of a test set point having the same class permutation as a training set point is at most $\frac{10^{6}}{10^{k}}$. For large $k$, we expect each test set point to test a class permutation unobserved in the training set.

Note that this task fits the modular structure of Equation 13: $y_{j}(x;U_{j})$ is a $1$-hot encoded label of dimensionality $10k$ indicating both the label of the $j$-th component image and which of the $k$ images is being predicted by the module. The full output is constructed as a sum of the $k$ labels $y_{j}(x;U_{j})$. As with the sine wave regression task, by increasing the number of component images, we may test generalization in arbitrarily high dimensions.

Figure: Figure 7: An example input from the Compositional CIFAR-10 task with $4$ component images. The goal of the task is to predict the classes of all component images in the input. In this case, with $4$ component images, the output target would be a $4$-hot encoded $40$-dimensional vector representing the true class of each of the $10$ possible classes for each component image.
Refer to caption: extracted/5841346/figs/compositional_cifar.png

#### Class combination experiments: Compositional CIFAR-10

We consider a Compositional CIFAR-10 variant in which the training inputs are constructed to have a distinct set of class label combinations compared to the test inputs (e.g. with $k=3$, if any training set input has the class combination cat, airplane, ship, then this class combination is not permitted on any test set input). This is done by partitioning the full set of class combinations into a set allocated for the training inputs and another disjoint set allocated for the test inputs. Thus, this tests out-of-distribution generalization. All other dataset parameters are set the same way as in the original Compositional CIFAR-10 task.

#### Noisy inputs experiment: Compositional CIFAR-10

We consider a Compositional CIFAR-10 variant in which the training inputs have added Gaussian noise drawn from $\mathcal{N}(0,\sigma I)$ of varying magnitude $\sigma$; this is done after the concatenation of $k$ images together. Test inputs do not have any added noise added. All other dataset settings are identical to the original Compositional CIFAR-10 task.

### E.2 Experiments on Monolithic Networks

#### Architecture and hyperparameter settings: sine wave regression

In our experiments, all neural networks are fully connected and use ReLU activations except at the final layer. We do not use additional operations in the network such as batch normalization. Networks are trained using Adam (Kingma & Ba, ) to minimize a mean squared error loss. We perform a sweep over learning rates in $\{0.001,0.01,0.1\}$ and find that the learning rate of $0.01$ performs best in general over all experiments, as justified by Fig 8 in Appendix E. Our results are reported in this setting. All networks are trained for $10000$ iterations which we find to be generally sufficient for convergence of the training loss.

The network architectures are varied as follows: the width of the hidden layers is selected from $\{8,32,128\}$, and the number of layers is selected from $\{3,5,7\}$. This yields $9$ total architectures. In order to consistently measure the number of parameters for an architecture as the input dimensionality varies, when we count the number of parameters we treat the input dimensionality $m$ as fixed at $m=1$. Note that this slightly underestimates the true number of parameters in each NN. The values of $k$ and $m$ range from $2$ to $9$. The value of $n$ ranges from $100$ to $100000$.

All experiments are run over $5$ random seeds and results are averaged. Experiments are run on a computing cluster with GPUs ranging in memory size from 11 GB to 80 GB.

#### Fitting our theoretical model: sine wave regression task

Our theoretical model of generalization has three free parameters $c$, $\Omega$ and $\alpha$. We select these parameters to best match the empirically observed trends of training set performance on the sine wave regression task; we find that $c=1.15$, $\Omega=1.57$, and $\alpha=0.85$.

#### Architecture and hyperparameter settings: Compositional CIFAR-10

In our experiments, all neural networks are fully connected and use ReLU activations except at the final layer. Note that the inputs are flattened; our networks do not use the spatial structure of the input. Batch normalization is applied before each ReLU. Images are normalized with standard normalization.

Networks are trained using Adam (Kingma & Ba, ) to minimize softmax cross-entropy loss using a learning rate of $0.0001$ and batch size of $128$. Note that the loss is averaged over all component images for each input. All networks are trained for a single epoch on $1000000$ random training points; note that with a moderately large number of images in each input, the total number of possible training inputs can be much larger. We use a test set of size $10000$.

The network architecture consists of fully connected layers of size $512$, $512$, $512$, $256$ and $128$ before a final fully connected layer to predict the output label.

Experiments are run over five random seeds for each hyperparameter configuration. Experiments are run on a computing cluster with GPUs ranging in memory size from 11 GB to 80 GB.

### E.3 Experiments on Modular Networks

#### Architecture and hyperparameter settings: sine wave regression task

Each module $\hat{y}_{j}(x;\hat{U}_{j})$ is constructed as follows: $\hat{U}_{j}$ consists of two vector components $\hat{u}_{j}$ and $\hat{v}_{j}$. The module output is constructed as:

$$ $\hat{y}_{j}(x;\hat{U}_{j})={f}_{j}(\hat{u}_{j}^{T}x)$ (72) $$

where $f_{j}$ represents a neural network with scalar input and output.

We set the kernel $\kappa$ as follows:

$$ $\displaystyle\kappa(x_{1},x_{2};(\hat{u},\hat{v}))$ $\displaystyle=e^{-\frac{1}{2\sigma^{2}}\left(\frac{x_{1}^{T}\hat{u}}{\hat{v}^{ T}\hat{u}}-\frac{x_{2}^{T}\hat{u}}{\hat{v}^{T}\hat{u}}\right)^{2}}$ (73) $\displaystyle+e^{-\frac{1}{2\sigma^{2}}\left(x_{1}-\frac{x_{1}^{T}\hat{u}}{ \hat{v}^{T}\hat{u}}\hat{v}-x_{2}+\frac{x_{2}^{T}\hat{u}}{\hat{v}^{T}\hat{u}} \right)^{2}\hat{v}}$ $$

where $\sigma$ is a hyperparameter. Intuitively, $\hat{v}\in\mathbb{R}^{m}$ corresponds to a direction along which projection directions $\hat{u}$ of other modules are not sensitive. This choice of kernel is motivated by the observation that if $\hat{v}^{T}\hat{u}_{j}$ for $j\neq i$, then $\left(x-\frac{x^{T}\hat{u}_{i}}{\hat{v}^{T}\hat{u}_{i}}\hat{v}\right)^{T}\hat{
u}_{j}=x^{T}\hat{u}_{j}$ and $\left(x-\frac{x^{T}\hat{u}_{i}}{\hat{v}^{T}\hat{u}_{i}}\hat{v}\right)^{T}\hat{
u}_{i}=0$: $x-\frac{x^{T}\hat{u}_{i}}{\hat{v}^{T}\hat{u}_{i}}\hat{v}$ removes all the information in $x$ relevant to module $u_{i}$ while retaining information relevant to all other modules.

Due to the computational cost of computing sample complexity via binary search (Algorithm 2), we fix several hyperparameters by small-scale experiments before the final experiment to control the total runtime. We first sweep through some parameters of in our module initialization method (Algorithm 1). Specifically, we perform sweep over module batch size in $\{32,128,512\}$, module learning rate in $\{0.001,0.01,0.1\}$, and module iteration number in $\{100,1000\}$. The combination of module batch size $128$, module learning rate $0.1$, and $1000$ module iterations has the lowest test error in our experiment.

Then, we sweep through number of architectural modules in $\{1,k,2\times k,5\times k\}$, learning rate in $\{0.001,0.01,0.1\}$ and $\sigma$ in $\{0.3,1.0,3.0\}$. We find that the combination of module number $5\times k$, learning rate $=0.01$ and $\sigma=1.0$ achieves best test performance. In addition, for our main experiments on modular NNs, all networks for dimension ($k=m$) smaller than 7 are trained for $1000$ iterations, while dimension 7 and 8 networks are trained for $700$ iterations, dimension 9 networks are trained for $500$ iterations and dimension 10 networks are trained for $200$ iterations. The high-dimension networks are stopped early since a smaller number of iterations was sufficient to converge on the training set; these numbers are determined based on small-scale experimental observations.

In binary search, we stop the search when the higher bound and the lower bound are close enough ($r-l<0.3$) to shorten our runtime. Also, due to GPU memory limitations, we can only support a sample size up to $10^{6}$, so we stop our experiments when the current sample size reaches $2^{22}$ ($c=22$). We set the maximum search iteration ($B$) to be $18$.

We test our network in 9 network architectures: the width of the hidden layers is selected from $\{8,32,128$}, and the number of layers is selected from $\{2,4,6$}. We also have three different values ($\{0.5,1.0,1.5$}) for the desired test error $\epsilon$ in binary search so as to pinpoint the most suitable value for further application of our method; we select a desired error of $1.5$ for our results. We keep $k=m$ in all experiments and the values range from $2$ to $9$. All experiments are run over $5$ random seeds and on the same computing cluster described in the previous section.

Figure: Algorithm 1 Finding a single module projection $\hat{U}$

Figure: Algorithm 2 Binary search for sample complexity

#### Disentanglement experiments: sine wave regression task

For our disentanglement experiments evaluating whether modular NNs can find the true modules underlying the task, we use the following hyperparameter settings: $m=k=10,n=1000$. We use a modular NN with $20$ modules. Each module uses a fully connected architecture with $5$ layers and a hidden layer width of $32$. We train the modular NN with a learning rate of $0.001$ for $1000$ iterations.

For learning our module initialization, we use $100$ iteration steps with a learning rate of $0.01$ and a batch size of $128$. $\sigma$ is set to $1.0$.

For constructing a t-SNE embedding, we use a perplexity of $5$. For computing similarity scores, we first compute the absolute value of the cosine similarity between each pair $(u_{i},\hat{u}_{j})$ of learned and target module directions. For each learned module direction $\hat{u}_{i}$, we then find the target module with the largest absolute cosine similarity. Finally, we average the maximum absolute cosine similarities across all modules to produce a similarity score: $\sum_{i=1}^{K}\max_{j=1}^{k}\frac{\hat{u}_{i}^{T}u_{j}}{||\hat{u}_{i}||_{2}||u
_{j}||_{2}}$.

#### Ablation experiments: sine wave regression task

For our ablation experiments, we train on $10000$ points. $k=m$ is varied from $2$ to $9$ and the number of architectural modules is set to $5\times k$. Each module uses a fully connected architecture with $5$ layers and a hidden layer width of $32$. We train the modular NN with a learning rate of $0.001$ for $1000$ iterations.

For learning our module initialization, we use $100$ iteration steps with a learning rate of $0.01$ and a batch size of $128$. $\sigma$ is set to $1.0$.

#### Architecture and hyperparameter settings: nonlinear sine wave regression task

To learn the nonlinear sine wave regression function, we consider modular architectures of the form: $\frac{1}{\sqrt{K}}\sum_{j=1}^{K}\hat{y}_{j}(||\hat{u}_{j}-x||_{2})$ where $\hat{y}_{j}$ is a neural network and $\hat{u}_{j}$ are learned parameters. We apply our method to learn an initialization for modules $\hat{u}_{j}$ using the following kernel:

$$ $\kappa(x_{1},x_{2};\hat{u}_{i})=e^{-\frac{1}{2\sigma^{2}}(||x_{1}-\hat{u}_{i}| |-||x_{2}-\hat{u}_{i}||)^{2}}$ (74) $$

We consider modules constructed as fully connected networks with $6$ layers and width $128$. We set $k=m=5$ and set $n=1000$. All other hyperparameter settings are consistent with our original sine wave regression experiments.

#### Architecture and hyperparameter settings: Compositional CIFAR-10

Note that given an input composed of $k$ images, the flattened input dimensionality is $3072k$. Each module $\hat{y}_{j}(x;\hat{U}_{j})$ is constructed as follows: $\hat{U}_{j}\in\mathbb{R}^{3072k\times 512}$ and the module output is constructed as:

$$ $\hat{y}_{j}(x;\hat{U}_{j})={f}_{j}(\hat{U}_{j}^{T}x)$ (75) $$

where $f_{j}$ represents a neural network with a $512$ dimensional input and an output of dimension $10k$.

We set the kernel $\kappa$ as follows:

$$ $\displaystyle\kappa(x_{1},x_{2};\hat{U})$ $\displaystyle=e^{-\frac{1}{2\sigma^{2}}\left|\left|{x_{1}^{T}\hat{U}}-{x_{2}^{ T}\hat{U}}\right|\right|_{2}^{2}}I$ (76) $$

where $\sigma$ is a hyperparameter; we set $\sigma=20.0$ to match the scale of the distances between projected inputs. To make kernel optimization more efficient, we stochastically optimize only the components of $y(X)^{T}\mathbf{K}^{-1}y(X)$ corresponding to a single class at a time.

Unless otherwise specified, hyperparameters are set to be consistent with the monolithic network. Additional hyperparameters are set as module batch size $128$ and module learning rate $0.01$. Module optimization is performed over a single pass over the training data. We fix the number of architectural modules as $32$. Each module is a ReLU-activated neural network with hidden layer sizes $256$, $128$ and $64$. Each ReLU is preceded by batch normalization. The module outputs are all concatenated and fed into a final linear layer to produce the $10k$ dimensional output. Critically, all the modules have the same weights. This is done to match the properties of the task: each modular component of the task is the same, namely to predict the class of a single CIFAR-10 image. This is unlike the sine wave regression task, where each modular component corresponds to a different function.

We vary the number of images $k$ from $1$ to $8$. All experiments are run on $5$ random seeds and on the same computing cluster described in the previous section.

## Appendix F Additional Experiments

Figure: Figure 8: Performance of monolithic architecture on a sine wave regression task across learning rates of $\{0.001,0.01,0.1\}$ (left), widths of $\{8,32,128\}$ (center), and number of layers $\{3,5,7\}$ (right). Each line represents a different architecture and desired test error, and the average performance is shown in bold. Note in the left plot that a learning rate of $0.01$ tends to perform the best, so all our experiments are reported in this setting.
Refer to caption: extracted/5841346/figs/monolithic_results_lr.png

Figure: Figure 9: Test loss of various training methods of training a modular NN on a sine wave regression task as the dimensionality of the task $k=m$ is varied. Margins indicate standard deviations over $5$ random seeds. Ours indicates that NN module directions $\hat{u}$ are initialized using our method and then trained on the task. Fixed initialization indicates that module directions are learned with our method and are fixed during task training. True modules indicates that $\hat{u}$ is set to the underlying module directions of the task $u$ (which are generally unknown) and fixed during task training. Observe that when $k=m$ is large, using the true, ground-truth module directions slightly outperforms our method, although interestingly our method performs better for low-dimensional tasks. Fixing the initialization found by our method results in significantly worse performance relative to allowing the module directions to vary.
Refer to caption: extracted/5841346/figs/ablations.png

Figure: Figure 10: Test loss of a modular NN initialized with our method with different numbers of modules on a sine wave regression task. The experiment is conducted with all hyperparameters fixed at the optimal value detailed in App E except the learning rate and $\sigma$. The experiment is run under 5 random seeds. The lines are averaged over all runs. Model performance increases as the number of modules increases. Intuitively, as the number of modules increases, the model can contain more information, thus achieving lower error.
Refer to caption: extracted/5841346/figs/n_modules.png

Figure: Figure 11: t-SNE embedding of target ($u$) and learned module projections ($\hat{u}$) learned by a randomly initialized modular NN without training (left), randomly initialized modular NN trained with gradient descent (center), and our initialized modules (right) on a sine wave regression task. From left to right, learned modules cluster more closely around target modules.
Refer to caption: x1.png

Figure: (a) Baseline Modular, Initialization
Refer to caption: x2.png

**Table 3: Comparison of our method with baselines on a Compositional CIFAR-10 variant in which the training inputs are constructed to have a distinct set of class label combinations compared to the test inputs (e.g., if any training set input has the class combination cat, airplane, ship, then this class combination is not permitted on any test set input). Thus, this tests combinatorial out-of-distribution generalization. We find test set accuracies for inputs with $6$ component images. Standard errors over $5$ trials are reported.**
| Method | Test Accuracy |
| --- | --- |
| Baseline monolithic | $42.56\%\pm 0.07\%$ |
| Baseline modular | $45.26\%\pm 0.06\%$ |
| Our method | $49.90\%\pm 0.15\%$ |

**Table 4: Comparison of our method with baselines on a Compositional CIFAR-10 variant in which the training inputs have added Gaussian noise drawn from $\mathcal{N}(0,\sigma I)$ of varying magnitude $\sigma$. We find test set accuracies for inputs with $6$ component images (note that the test points do not have added noise). This tests out-of-distribution generalization to small distribution shifts. Standard errors over $4$ trials are reported.**
| Noise Level | Baseline Monolithic | Baseline Modular | Our Method |
| --- | --- | --- | --- |
| 0 | $42.92\pm 0.05\%$ | $45.66\pm 0.14\%$ | $50.49\pm 0.09\%$ |
| 0.3 | $42.57\pm 0.05\%$ | $45.42\pm 0.10\%$ | $50.27\pm 0.08\%$ |
| 1 | $39.47\pm 0.09\%$ | $43.07\pm 0.08\%$ | $46.95\pm 0.11\%$ |
| 3 | $31.75\pm 0.03\%$ | $34.71\pm 0.05\%$ | $34.67\pm 0.13\%$ |
| 10 | $21.76\pm 0.10\%$ | $24.31\pm 0.34\%$ | $24.41\pm 0.13\%$ |
| 30 | $11.09\pm 0.33\%$ | $12.50\pm 0.55\%$ | $12.29\pm 0.40\%$ |

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 3
presentation: 2
questions: - What is the explicit definition of modularity being used? This concept was referenced without ever being really explicitly defined in a general-but-still-technical sense (and attention was given as an example). I ended up being confused about whether the focus was on independently functioning parts of a network, or shared weights in a more general sense .
- As mentioned in "Weaknesses": what assumptions are made in general about the structure of the feature matrix? It is indicated in the modular version of the model that arbitrary nonlinear transforms of the input are considered as valid feature matrices, but this isn't clarified for the first treatment of the model 
- Do the benefits cited require being correct about the number (k) of modules in the underlying data? How much do positive results depend on being correct in your choice of k relative to what is present in the underlying data?
soundness: 3
strengths: - The paper focuses on an empirical task that lets them validate their ideas, but is realistic in terms of scope 
- The paper engages with an interesting problem of the structure of weights that allow for better generalization
summary: This paper (1) constructs a simplified theoretical model of generalization (focused on the case of linear regression from what I believe to be a set of not-strictly-defined features), (2) provides empirical demonstrations that, at least in broad strokes, a sine wave regression task follows the predictions made by that model, (3) argues that in cases where the simple problem under consideration is modular (here meaning that it has k modules of size b which interact, rather than a full set of P features), better sample complexity can be obtained by using an explicitly modular parameter structure, and (4) provides some empirical support for this generalization behavior on another set of (this time modular) sine wave regression tasks.
weaknesses: - The paper's notation, particularly in the initial presentation of the theoretical model, was confusing and felt under-explained. In particular I was confused by feature matrix (how did it get constructed from the inputs? What assumptions are being made about it? Why is it a matrix to begin with rather than a feature vector), and this confusion made it hard to understand future claims made in the paper (especially since the central claim was about the effect on generalization of _input_ dimension, which is mediated by the function implied in this feature matrix)
- The forms of the expected training and test loss could have been broken down in a clearer and more intuitive way, rather than simply being presented as not-very-comprehensible formulas 
- This paper assumes that the only way to benefit from the generalization behavior of modularity is to have explicitly modular structure; it would have been interesting if it had also engaged with whether modular data gives you generalization benefits without a parallel parameter structure (since in practice modern models seem to generalize well without the benefit of this)

---

confidence: 2
contribution: 3
initial_rating: 8
presentation: 4
questions: Continuing with the previous weakness, the algorithm to learn or initialize the input projection parameters may not work in such a case as it is dependent on the initial module NN weights. 

The generalization performance for the compositional CIFAR-10 experiment can be divided into input class permutations present in the training data vs. not present in the training data to further dissect the difference between monolithic NNs and modular NNs. The sample complexity experiments with compositional CIFAR-10 tasks are not present and should be added to further strengthen the claims. 

For individual tasks considered, there appears to be a large amount of tuning of methodology to train the modular NNs and the module input projections. (Referring to appendix)

How would the solution to training and test errors change if the input projections from each module were removed, and the modules considered the input x in its entirety? This is consistent with the current mixture-of-exert (MoE) models. 
 
Is there a validation set used for the experiments or is the generalization performance reported from the last training iteration?

Do the modular NNs treat the number of modules as a hyper-parameter and tune it to improve the performance. Or is it an architectural characteristic such as the width or depth in monolithic NNs that is fixed ? 

The CIFAR-10 experiments are run only for a single epoch, will increasing the number of epochs result in better performance for networks? 


Minor: Equations are referred to the appendix when they are also present in the main part of the paper.
soundness: 4
strengths: The theoretical derivations are sound and very well done, and the paper is well written. The authors did a good job showcasing how modular networks can be related to monolithic networks (based on certain linear assumptions).

Intuitively, modeling of the data, and learning of parameters related to module input bottlenecks make sense. This is similar to learning the connectivity associated with module input or limiting the number of module inputs to avoid module collapse. 

The experiments show a clear trend that the input projection mechanism results in better performance and sample complexity, as compared to monolithic NNs and end-to-end trained modular NNs.
summary: The paper shows that the sample complexity associated with training modular neural networks is independent (under certain conditions) of the input dimensionality and does not follow the same exponential increase with input dimension as in the case of monolithic or traditional neural networks (NNs). 

First, the authors present a derivation of training and test error in monolithic NNs when the task is non-modular. The task and the NN are modeled linearly based on features generated using the input. The authors then, empirically validate that NNs with different architectures (loosely) match the theoretical trends when varying the input dimension, number of samples and the number of parameters. Note that the task considered for this experiment is modular. 

The authors then theoretically compute the training and generalization errors of modular NNs given that the underlying task is also modular with the same structure. Each module is associated with a small NN (modeled linearly) and an input projection mechanism that reduces the module input dimensionality. The task is also modeled in the same way where the parameters are randomly initialized. 

The resulting closed form solution shows that the training error is independent of the input dimension and the test error under the condition of under-parametrization is also independent of the input dimension. This result hinges on the input projection associated with each module, where the dimensionality associated with the overall input is reduced. The authors then propose a method to initialize (or learn) the parameters associated with the input projections. Once initialized the modular NNs are trained end-to-end. 

Empirical results show that modular NNs that are learned using the proposed initialization mechanism achieve significant improvements over monolithic NNs and modular NNs conventionally trained.
weaknesses: The major weakness of the paper is the consideration of a single layer of modules and data generating system. In such a system the output is a linear composition of the module outputs. This may not be true for many real world systems where multiple such modular layers can exist in a hierarchy.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: (1) I was confused by the results in figure 1, whereby test error INCREASES with larger datasets. This seems incompatible with empirical observartions and traditional statistical analyses of generalization.

(2) I did not understand in what sense the y_j could be considered independent (and what is the random variable), after eq 3.

(3) eqn 4 seems wrong: on the LHS y is a function of the linear projection U x, whereas on the RHS U and x only interact via the presumably non-linear function phi.

(4) why should we expect eqn 15 to give the minimum norm solution? (i.e. why is it a solution and why is it minimum norm, from what class)

(5) You should add the citation of the original mixture of expert papers (e.g. Jacobs et al 1991).
soundness: 2
strengths: Much remains to be understood about the generalization behavior of neural nets, especially the types that have a modular architecture, so advances on the theory (in special cases) that are presented do seem useful.
summary: The authors present both theoretical and algorithmic results regarding the generalization ability of a particular form of modular architectures. They first highlight theoretical results on generalization results showing exponentially large sample complexity as a function of the input dimension. Then they present a particular class of modular architectures as a sum of experts and assume that the training data was generated by this same architecture. They show that thanks to each module making a low-dimensional projection before processing the data, the scaling behavior is better behaved. Then they present a kernel-style algorithm to initialize such an architecture, to be then fine-tuned by usual supervised learning and SGD. Finally, they show results on a toy 1-dim sine-wave regression task and on the recently introduced compositional MNIST task.
weaknesses: (1) the theoretical results are not suprising: projecting the m-dim input to several b-dim low-dimensional representations unsurprisingly reduces the exponential badness from m to b. The theory is also of fairly limited scope, with lots of unreasonable assumptions (e.g., of linearity wrt parameters) that may not tell us as much as we would like for more general forms of modular architectures.

(2) the proposed algorithm is unlikely to scale well in terms of computational efficiency beyond small-size problems and into frontier AI, given the use of kernel methods in the novel part of the method

(3) the fact that all modules are initialized independently and using the same (randomized) procedure suggests that a significant part of the advantage could come from an ensemble effect (which always helps generalization)

(4) the paper seems to overclaim in multiple places, e.g., suggesting that their results tracks empirical behavior of modern neural nets (even the empirical comparisons don't match the theory, e.g., fig 2 bottom right).

(5) I did not find numerical comparisons against benchmark results from other papers, and when I look at Jarvis et al 2023, their figures show much lower errors. Hence the experimental results may not be that good after all.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: See above
soundness: 4
strengths: This paper presents an excellent set of theoretical results outlining expected scaling laws for modular networks when task have a modular structure. It also presents a practical initialization scheme for potential models that is based on self-supervised alignment with task statistics. Lastly, it validates predictions with non-trivial and relevant experiments with architectures that go beyond the ones used for theory, outlining the potential generality of the result.
summary: In this paper, the authors study modular neural networks and their ability to learn functions that are themselves modular insofar as they have either a compositional nature or a modular composition in their statistics. The authors first provide a theorem outlining expected scaling laws for sample complaxitiy and trainign accuracy in a curated linear task setup. They then corroborate the predictions from this theorem with numerical experiments. Following this, they present a task setup with explicit ground-truth modular structure and prove a second theorem outlining similar sclaing propoerties for this novel tasks setup along with a particual modular NN architecture. The result is reduced scaling complexity for tusch tasks when NNs are appropriately modular. The authors then propose an initialization scheme for NNs which first learn module initialization in a self-supervised fashion using task statistics. Finally, the authors show that this approach behaves as expected on non-trivial tasks such as compositional CIFAR, and that the proposed method works well on other modular architectures beyond the one used for their theoretical results.
weaknesses: The paper is sound. There are potntials for improvements in two key areas:

1. In experiments, it is unclear if the generalization advanatage of the modular networks remain if one factors in the pre-training (i.e. learning module initialization). In other words, for the same total compute, would a monolithic model do as well as the modular one for which some of the compute budget went toward initialization? My apologies if I missed it, but this is a key point that would factor into real-world scaling laws.

2. What happens if there is a missmatch between the task modularity and the NNm odule count? In experiments, the models could accurately recover tasks modularity when the number of modules is known. In contrast, in other tasks, this is not known a priori. How would the models behave in this mismatched environment? To be fair, the authors acknowledge this issue, but I wonder if some rapid experiments could outline expected behaviors in this case. Once again, apologies if this has been explored and I missed it.

3. While this is no doubt very relevant for the field, it's overall impact with respect to modern architectures remains unclear. Some discussion about the use of such methods in modern settings such as in the prsence of attention could greatly enhance the scope of the result. This is not necessary however, as this is a good paper in the current scope.

Meta-review: This paper derives a theoretical quantitative analysis of the generalization abilities of modular architectures. The analysis shows that under certain assumptions (linearity, specific modular structure of architecture and data) models can avoid exponential sample complexity associated with high dimensional inputs. The paper corroborates this with empirical result on two toy tasks.

Reviews agree that the studied problem is interesting, but disagree substantially about the relevance of the results. The main concern is that the assumptions might be too restrictive to apply to realistic cases and modern architectures. The authors defend their paper as a stepping stone towards better understanding of modular generalization. 

From the reviews it seems clear that the results will be of interest to at least a part of the community. I lean towards counting being controversial as a potential advantage in this case, and therefore recommend accepting this paper.

Decision: Accept (Poster)