The Persian Rug: solving toy models of superposition using large-scale symmetries

## Abstract

Abstract We present a complete mechanistic description of the algorithm learned by a minimal non-linear sparse data autoencoder in the limit of large input dimension. The model, originally presented in Elhage et al. ( 2022 ) , compresses sparse data vectors through a linear layer and decompresses using another linear layer followed by a ReLU activation. We notice that when the data is permutation symmetric (no input feature is privileged) large models reliably learn an algorithm that is sensitive to individual weights only through their large-scale statistics. For these models, the loss function becomes analytically tractable. Using this understanding, we give the explicit scalings of the loss at high sparsity, and show that the model is near-optimal among recently proposed architectures. In particular, changing or adding to the activation function any elementwise or filtering operation can at best improve the model’s performance by a constant factor. Finally, we forward-engineer a model with the requisite symmetries and show that its loss precisely matches that of the trained models. Unlike the trained model weights, the low randomness in the artificial weights results in miraculous fractal structures resembling a Persian rug, to which the algorithm is oblivious. Our work contributes to neural network interpretability by introducing techniques for understanding the structure of autoencoders. Code to reproduce our results can be found at https://github.com/KfirD/PersianRug .

## 1 Introduction

Large language model capabilities and applications have recently proliferated. As these systems advance and are given more control over basic societal functions, it becomes imperative to ensure their reliability with absolute certainty. Mechanistic interpretability aims to achieve this by obtaining a concrete weight-level understanding of the algorithms learned and employed by these models. A major impediment to this program has been the difficulty of interpreting intermediate activations. This is due to the phenomena of superposition, in which a model takes advantage of sparsity in the input data to reuse the same neurons for multiple distinct features, obscuring their function. Finding a systematic method to undo superposition and extract the fundamental features encoded by the network is a large and ongoing area of research Bricken et al. (); Cunningham et al. (); Gao et al. (); Engels et al. ().

Currently, the most popular method of dealing with superposition is dictionary learning with sparse autoencoders. In this method, the smaller space of neuron activations at a layer of interest is mapped to a larger feature space. The map is trained to encourage sparsity and often consists of an affine + ReLU network. This method has been applied to large language models revealing many strikingly interpretable features (e.g. corresponding to specific bugs in code, the golden gate bridge, and sycophancy), even allowing for a causal understanding of the model’s reasoning in certain scenarios Marks et al. ().

The sparse decoding ability of the affine + ReLU map was recently studied in the foundational work Elhage et al. (), which introduced and studied a toy model of superposition. The model consisted of a compressing linear layer modeling the superposition(^1^11This is because, if good enough recovery is possible for most features, the pigeonhole principle tells us that at least some of the smaller space activations must encode information about multiple input features.) followed by a decompressing affine + ReLU layer, trained together to auto-encode sparse data. They showed that the network performs superposition by encoding individual feature vectors into nearly orthogonal vectors in the smaller space. The affine layer alone is unable to decode sparse input vectors sufficiently well to make use of superposition, but the addition of the ReLU makes it possible by screening out negative interference.

While Elhage et al. () provides valuable empirical and theoretical insights into superposition, it does not obtain a strong enough description of the model algorithm to quantitatively characterize the algorithm’s performance. Given the extensive use of the affine + ReLU map for decoding sparse data in practice, it is important to obtain a complete analytic understanding of the model behavior over a large parameter regime. As we will see, this will inform the design of better sparse autoencoder architectures.

In this work we obtain such an understanding by considering a particularly tractable regime of the Elhage et al. () model: permutation symmetric data (no input feature is privileged in any way), and the thermodynamic limit (a large number of input features), while maintaining the full range of sparsity and compression ratio values. In this regime, the learned model weights are permutation symmetric on large scales, which sufficiently simplifies the form of the loss function to the point where it is analytically tractable, leaving only a small number of free parameters. We then forwards-engineer an artificial set of weights satisfying these symmetries and optimizing the remaining parameters, which achieves the same loss as a corresponding trained model, implying that trained models also implement the optimal permutation symmetric algorithm. The artificial set of weights resembles a Persian rug fig. 1, whose structure is a relic of the minimal randomness used in the construction, illustrating that the algorithm relies entirely on large-scale statistics that are insensitive to this structure. Finally, we derive the exact power-law scaling of the loss in the high-sparsity regime.

We expect our work to impact the field of neural network interpretability in multiple ways. First, our work provides a basic theoretical framework that we believe can be extended to other regimes of interest, such as structured correlations in input features, which may help predict scaling laws in the loss based on the data’s correlations. Second, our work rules out a large class of performance improvement proposals for sparse autoencoders. Finally, our work provides an explicit example of a learned algorithm that is insensitive to microscopic structure in weights, which may be useful for knowing when not to analyze individual weights.

The paper is structured as follows. In section 2 we review the model and explain our training procedure. In section 3 we show empirically that large models display a “statistical” permutation symmetry. In section 4 we extract the algorithm by plugging the symmetry back into the loss, introduce the Persian rug model which optimizes the remaining parameters, show that large trained models achieve the same loss, and derive the loss behavior in the high sparsity limit. In section 5.2 we discuss related works, and conclude with section 5.

## 2 The Model

We consider the following model for our non-linear autoencoder with matrix parameters $W_{\text{in}}\in\mathbb{R}^{n_{d}\times n_{s}},W_{\text{out}}\in\mathbb{R}^{n_
{s}\times n_{d}},{\mathbf{b}}\in\mathbb{R}^{n_{s}}$,

$$ $\displaystyle f_{\text{nonlinear}}({\mathbf{x}})$ $\displaystyle=\operatorname{ReLU}(W_{\text{out}}W_{\text{in}}{\mathbf{x}}+{ \mathbf{b}}).$ (1) $$

$W_{\text{in}}$ is an encoding matrix which converts a sparse activation vector ${\mathbf{x}}$ to a dense vector, while $W_{\text{out}}$ perform the linear step of decoding. We also consider a simple model for the sparse data on which this autoencoder operates. Each vector ${\mathbf{x}}$ is drawn i.i.d. during training, and each component is drawn i.i.d. with $x_{i}=c_{i}u_{i}$, where $c_{i}\sim\text{Bernoulli}(p)$ and $u_{i}\sim\text{Uniform}[0,1]$ are independent variables. This ensures that ${\mathbf{x}}$ is sparse with typically only $pn_{s}$ features turned on.

We train our toy models to minimize the expected $L_{2}$ reconstruction loss,

$$ $\displaystyle L({\mathbf{x}};W_{\text{out}},W_{\text{in}},{\mathbf{b}})=E||{ \mathbf{x}}-f_{\text{nonlinear}}({\mathbf{x}})||_{2}^{2}.$ (2) $$

It is known that for the linear model (eq. 1 without the ReLU), the optimal solution is closely related to principle component analysis (see, for example, Plaut () and p. 563 or Bishop & Nasrabadi ()). In particular, the reconstruction loss decreases linearly in the hidden dimension $n_{d}$ when all features are i.i.d. On the other hand, the model eq. 1 will have a much quicker reduction in loss, as will be described in section 3.1.

For all models presented we train them with a batch size of $1024$ and the Adam optimizer to completion. That is training continues as long as the average loss over the past 100 batches is lower than the average loss over the 100 batches prior to that one. Our goal with training is to ensure that we have found an optimal model in the large-data limit to analyze the structure of the model itself.

## 3 Empirical observations

In this section, we present empirical observations of the trained models. We start by presenting a remarkable phenomenon this model exhibits in the high-sparsity regime: a dramatic decrease in loss as a function of the compression ratio. We then turn to a mechanistic interpretation of the weights which gives empirical evidence for the phenomena needed to understand the algorithm the model learns. These are manifestations of a partially preserved permutation symmetry of the sparse degrees of freedom.

### 3.1 Fast loss drop

To gauge the performance of the model, we plot the loss (eq. 2) as a function of the compression ratio $\mathchoice{\frac{n_{d}}{n_{s}}}{n_{d}/n_{s}}{n_{d}/n_{s}}{n_{d}/n_{s}}$. In fig. 2 we plot the performance of the linear versus non-linear models for representative parameters. It is clear that the non-linear model outperforms the linear model up until near the $\mathchoice{\frac{n_{d}}{n_{s}}}{n_{d}/n_{s}}{n_{d}/n_{s}}{n_{d}/n_{s}}\approx
1$ regime due to an immediate drop in the loss.The slope and duration of this initial fall is controlled by $p$. In particular, in the high-sparsity regime ($p$ close to zero), the loss drops to zero entirely near the $\mathchoice{\frac{n_{d}}{n_{s}}}{n_{d}/n_{s}}{n_{d}/n_{s}}{n_{d}/n_{s}}\approx
0$ regime. What is going on here? To explain this behavior, we analyze the algorithm the model encodes.

Figure: Figure 2: Loss curves of trained models, Persian rug models, and optimal linear models as a function of the compression ratio.
Refer to caption: extracted/5946829/images/line_plot_target_p_feat_0.04.png

### 3.2 Statistical permutation symmetry

Rather than looking individually at the weights, it is helpful to look at the matrix $W=W_{\text{out}}W_{\text{in}}$, shown in fig. 3. As will be made precise in section 4, we are interested in the structure of each row since it is responsible for the corresponding feature output. Though the matrix is not permutation symmetric in detail, its statistics are sufficiently so as far as the algorithm is concerned. In particular, in the $\mathchoice{\frac{n_{d}}{n_{s}}}{n_{d}/n_{s}}{n_{d}/n_{s}}{n_{d}/n_{s}}\rightarrow\infty$ limit, the matrix is statistically permutation invariant in the following sense:

- 1.
the diagonal elements become the same (fig. 4),
- 2.
the bias elements become uniformly negative, which can be seen in the uniformity and slight blue shade in fig. 3 and is quantified in fig. 5,
- 3.
the statistics of the off-diagonal terms in each row are sufficiently uniform for interference to become Gaussian (fig. 7), and finally
- 4.
the mean and variance of each of these becomes the same across rows (fig. 6).

Figure: Figure 3: Plot of the first $30\times 30$ $W$ elements and the corresponding bias (${\mathbf{b}}$) components, at $p=4.5\%$ and ratio $n_{s}=512$. The diagonal components are all at similar values of $1.29\pm.01$ (one standard deviation) while the off-diagonal components are approximately mean-zero, appearing like noise. The bias elements are all negative around $-.18\pm.01$. This statistical uniformity is a permutation symmetry across the sparse features.
Refer to caption: extracted/5946829/images/model_weights_non_average_run_09_26_921.png

Immediately on looking at the exemplar model from fig. 3 we see that the diagonals are large and uniform while the off-diagonal components are fluctuating and much smaller. This serves the function of routing each input to the corresponding output, except with some off-diagonal contributions which should be as small as possible.

Figure: Figure 4: Permutation symmetry of diagonal values. We plot the mean-square fluctuation of the diagonal values corresponding to each model. Models are trained as a function of $p$ and $\mathchoice{\frac{n_{d}}{n_{s}}}{n_{d}/n_{s}}{n_{d}/n_{s}}{n_{d}/n_{s}}$. The emergence of symmetry as $n_{s}$ grows (at all locations in the diagram) is a crucial element of the algorithm implemented by the autoencoders.
Refer to caption: extracted/5946829/images/RPS_diag_var.png

Because the diagonal components appear uniform, we quantify this by looking at the fluctuation on a per-model basis. For each model, we compute the empirical mean and variance of the diagonal components and plot the latter in fig. 4. Initially, when $n_{s}=128$ is small, we see a region with large fluctuations among the diagonal components, particularly when the $p$ and $\mathchoice{\frac{n_{d}}{n_{s}}}{n_{d}/n_{s}}{n_{d}/n_{s}}{n_{d}/n_{s}}$ are both small. We hypothesize this occurs because for small models (and small feature probabilities) the contribution to any output from off-diagonal terms may fluctuate because of the small size of the matrix.

Figure: Figure 5: Permutation symmetry of bias values. We plot the mean-square fluctuation of values in the bias vectors corresponding to each model, which are trained as a function of $p$ and $\mathchoice{\frac{n_{d}}{n_{s}}}{n_{d}/n_{s}}{n_{d}/n_{s}}{n_{d}/n_{s}}$. As $n_{s}$ increases the fluctuation over bias elements generally decreases in all trained models.
Refer to caption: extracted/5946829/images/RPS_bias_var.png

We see similar features among the entries of the bias in fig. 5. The relatively large initial fluctuations generally decay as $n_{s}$ grows larger. Similar features are present in the two plots, such as the relatively low fluctuation regions near $(p,\mathchoice{\frac{n_{d}}{n_{s}}}{n_{d}/n_{s}}{n_{d}/n_{s}}{n_{d}/n_{s}})
\approx(.1,.5)$ because the bias and the diagonal entries are coupled when a diagonal element in a row is large the corresponding bias must be scaled up as well and vice versa.

For each of these algorithm-revealing features, we give three plots, each with increasing $n_{s}$, to show that the feature becomes apparent in the thermodynamic limit. Each plot is a heat map with the variable representing the feature plotted in color and spatial dimensions covering the entire parameter space of $\mathchoice{\frac{n_{d}}{n_{s}}}{n_{d}/n_{s}}{n_{d}/n_{s}}{n_{d}/n_{s}}$ and $p$.

Figure: Figure 6: Permutation symmetry of off-diagonal statistics. The symmetry breaking parameter $\Delta\operatorname{var}(\nu)$ is given by the variance across all rows of the squared sum of the off diagonal elements in each row, up to a constant. Once $n_{s}$ reaches $8192$ all noise variables have nearly identical variances.
Refer to caption: extracted/5946829/images/RPS_chi_varvar.png

We similarly see in fig. 6 that the off-diagonal terms concentrate as well. To be more precise, we look at the fluctuation of the total norm of the off-diagonal elements in each row in $W$ across rows. The norm is a measure of the amount of signal arriving at a particular output from other inputs. A smaller contribution from other inputs means that each output has a higher signal-to-noise ratio as long as the diagonal term remains constant.

Figure: Figure 7: Interference Gaussianity in the large $n_{s}$ limit. We show that the Lyapunov condition is better satisfied as $n_{s}$ grows larger, by plotting the worst-case Lyapunov condition (see eq. 3) over all rows of $W$. Models are trained at different $p$ and $\mathchoice{\frac{n_{d}}{n_{s}}}{n_{d}/n_{s}}{n_{d}/n_{s}}{n_{d}/n_{s}}$ values. We see that for small $n_{s}$ only models trained at large $p$ have nearly-uniform off-diagonal entries whereas all models approach uniformity at large $n_{s}$.
Refer to caption: extracted/5946829/images/RPS_Max_Lyapunov.png

Given that there are many off-diagonal terms, and those terms are all added together in the process of multiplying $W{\mathbf{x}}$, it is natural to ask if they satisfy the Lyapunov condition, which would imply that the contribution to output $i$ from off-diagonal terms $j\neq i$ is Gaussian. In fig. 7 we plot the statistic

$$ $\Lambda\equiv\max_{i}\frac{\sum_{j\neq i}|W_{ij}|^{3}}{\left(\sum_{j\neq i}W_{ ij}^{2}\right)^{3/2}}$ (3) $$

where the sum is over all values of $j\neq i$ for each $i$. The statistic tends towards 0 as $n_{s}$ grows larger which means that the contribution from the off-diagonal elements tends towards a Gaussian distribution. This is strong numerical evidence that no small set of terms in the off-diagonals dominate.

### 3.3 Optimization of residual parameters

The statistical permutation symmetry places constraints on the possible values of $W$ and ${\mathbf{b}}$. The constraint on ${\mathbf{b}}$ is straightforward: it is proportional to the all ones vector, i.e. there is a number $b$ such that ${\mathbf{b}}_{i}=b$ for all $i$. As will be explained in section 4.1, the relevant degrees of freedom remaining in $W$ are one number $a$ equal to the diagonals, and $\sigma$ characterizing the root mean square of the off diagonal rows. The precise values of the off diagonals can be thought of as irrelevant “microscopic information”. Thus there are three relevant degrees of freedom remaining: $b,a$, and $\sigma$.

In section 4.2.2 we give a specific set of values for $W_{\text{in}}$ and $W_{\text{out}}$ via the “Persian rug” matrix, which have the statistical permutation symmetry in the $n_{s}\rightarrow\infty$ limit while also optimizing $\sigma$. The remaining parameters, $a$ and $b$ can be optimized numerically. In fig. 2 we compare the loss curve of this artificial model with that of a trained model, and see that they are essentially the same.

## 4 Extracting the Algorithm

In this section we give a precise explanation of the algorithm the model performs. We start with a qualitative description of why the statistical permutation symmetry gives a good auto-encoding algorithm when the remaining macroscopic degrees of freedom are optimized. We then find an artificial set of symmetric weights with optimized macroscopic parameters. We show that the trained models achieve the same performance as the artificial model, thus showing they are optimal even restricting to statistically symmetric solutions. Finally, we derive an explicit form of the loss scaling and argue that ReLU performs near optimally among element-wise or “selection” decoders.

### 4.1 Qualitative Description

A key simplification is to consider strategies as collections of low-rank affine maps rather than as the collection of weights directly. In other words consider the tuple $(W,{\mathbf{b}})$ where $W=W_{\text{out}}W_{\text{in}}$ to define the strategy. We must restrict to $W$ with rank no more than $n_{d}$ because it is the product of two low-rank matrices. Given any such $W$ we may also find $W_{\text{in}}$ and $W_{\text{out}}$ of the appropriate shape (e.g. by finding the SVD), so the two representations are equivalent.

Let us now take a close look at the reconstruction error. Let us consider output $i$, which is given by

$$ $(f_{\text{nonlinear}}({\mathbf{x}}))_{i}=\operatorname{ReLU}(W_{ii}x_{i}+\sum_ {j=1,j\neq i}^{n_{s}}W_{ij}x_{j}+{\mathbf{b}}_{i})=\operatorname{ReLU}\left(a( x_{i}+\nu_{i})+{\mathbf{b}}_{i}\right)$ (4) $$

where we used that all diagonal elements of $W$ are equal to $a$ and defining $\nu_{i}=a^{-1}\sum_{j=1,j\neq i}^{n_{s}}W_{ij}x_{j}$. At this stage we can already see that the individual off-diagonal elements are only important insofar as they characterize the Gaussian random variable $\nu_{i}$. More specifically, because $\nu_{i}$ is completely characterized by its mean and variance, only the sum and sum of squares of the off-diagonal components matter. Because we have a bias term the mean of $\nu_{i}$ may be absorbed into it, so from now we assume that $\nu_{i}$ are all zero-mean. This simplification is why permutation symmetry emerges when $n_{s}$ is large, and allows us to consider only the macroscopic variables introduced when we discuss the loss.

The reconstruction error is (dropping the $i$ subscript, and substituting ${\mathbf{b}}_{i}=b$)

$$ $\displaystyle L=\mathbb{E}_{x,\nu}[(x-\operatorname{ReLU}\left(a(x+\nu)+b \right))^{2}].$ (5) $$

We can further decompose this by taking the expectation value over whether $x$ is on or off, so

$$ $\displaystyle L=(1-p)L_{\text{off}}+pL_{\text{on}}$ $$

where

$$ $\displaystyle L_{\text{off}}$ $\displaystyle=\mathbb{E}_{\nu}[\left(\operatorname{ReLU}\left(a\nu+b\right) \right)^{2}],\text{ and}$ $\displaystyle L_{\text{on}}$ $\displaystyle=\mathbb{E}_{u,\nu}[(u-\operatorname{ReLU}\left(a(u+\nu)+b\right) )^{2}]$ $$

with ${u\sim\text{Uniform}[0,1]}$. Let us explore the regimes of macroscopic parameters $a,b,\sigma$ when either or both of these loss terms are low.

The impact of all the non-diagonal terms has been summed up in the “noise” $\nu$. Though it is a deterministic function of ${\mathbf{x}}$, output $i$ has no way to remove $\nu$ from $u+\nu$. The best it can do is to estimate $u$ from $u+\nu$ because it doesn’t have any other information. This exemplifies a key principle – by restricting the computational capacity of our decoder, deterministic, but complicated correlations act like noise.

The main advantage of the nonlinear autoencoder is that the dominant contribution to the loss, $L_{\text{off}}$, can be immediately screened away by making $a\nu+b$ either small or negative, allowing the network to focus on encoding active signals. This immediate screening is always possible either by choosing $b+\mu$ large and negative or $a$ and $\mu+b$ small. However, these strategies come at a cost because the output value is distorted from $u$ to $au+b$. It is thus preferable instead for $\nu$ and $b$ to be as small as possible, which occurs when $\sigma$ is as small as possible. As we will see, $\sigma$ is the only parameter we are not free to choose in the large $n_{s}$ limit, and whose value will be bounded as a function of $\mathchoice{\frac{n_{d}}{n_{s}}}{n_{d}/n_{s}}{n_{d}/n_{s}}{n_{d}/n_{s}}$ and $p$. Since $L_{\text{off}}$ is the dominant contribution to the loss, therefore, it will thus be necessary to damp the signal by setting $a$ small and/or $b$ large and negative in regimes where $\sigma$ is uncontrollably large.

Given that we see a statistical permutation symmetry in trained models let’s consider symmetric strategies so that $W_{ii}=a$ and ${\mathbf{b}}_{i}=b$ for all $i=1,\dots,n_{s}$. We will show that optimizing the remaining macroscopic parameters makes $f_{\text{nonlinear}}$ act close to the identity on sparse inputs.

### 4.2 Optimizing the macroscopic parameters

We have seen qualitatively that a statistically symmetric strategy exists in certain regimes of the macroscopic parameters. Two of these parameters, $a$ and $b$ are unconstrained. Furthermore the loss should be monotonically increasing with $\sigma$ because a larger $\sigma$ implies more noise which hinders reconstruction. Thus we now prove lower bounds on $\sigma$ and construct an artificial set of statistically permutation symmetric weights which achieve this bound. Finally we will compare the reconstruction loss of this strategy with the learned one to justify that those ignored microscopic degrees of freedom were indeed irrelevant.

#### 4.2.1 Optimal σ 𝜎 \sigma italic_σ

Assuming the permutation symmetry we discovered earlier in our empirical investigations we will derive a bound on the variance of the output. Additionally for an optimal choice of ${\mathbf{b}}$ the average loss is increasing in the variance, because a larger variance corresponds to a smaller signal-to-noise ratio. Taken together these two facts will give a lower bound on the loss. We will then provide an explicit construction which achieves this lower bound and illustrates how the algorithm works.

The lower bound on the variance comes from the fact that $W$ is low-rank with constant diagonals. For now let us ignore the overall scale of $W$, and just rescale so that the diagonals are exactly 1. The bound we are about to prove is very similar to the Welch bound (Welch, ) who phrased it instead as a bound on the correlations between sets of vectors. We produce an argument for our context, which deals with potentially non-symmetric matrices $W$, the details of which are located in appendix C.

We show that

$$ $\sigma^{2}\geq\frac{4p-3p^{2}}{12}\left(\frac{n_{s}}{n_{d}}-1\right)$ (6) $$

with equality only when $W$ is symmetric, maximum rank, with all non-zero eigenvalues equal. This naturally leads to a candidate for the optimal choice of $W$, namely matrices of the form

$$ $W\propto OPO^{T}\text{ and }W_{ii}=1$ (7) $$

where $O$ is an orthogonal matrix and $P$ is any rank-$n_{d}$ projection matrix. This kind of matrix saturates the bound because it is symmetric and has all nonzero eigenvalues equal to 1.

A note on the connections between these matrices and tight frames – if we take the “square root” of $W$ as a $n_{s}\times n_{d}$ matrix $\sqrt{W}$ such that $\sqrt{W}\sqrt{W}^{\dagger}=W$ then the $n_{d}$ dimensional rows of $\sqrt{W}$ are a tight-frame. This is because $\operatorname{Tr}(W)^{2}=n_{s}^{2}=n_{d}\operatorname{Tr}WW^{\dagger}$ which is the variational characterization of tight-frames as in Theorem 6.1 from Waldron ().

#### 4.2.2 Persian rug model

We now give an explicit construction of an optimal choice for $W$. The construction is based on the Hadamard matrix of size $n_{s}=2^{m}$ for some integer $m$, defined as

$$ $\displaystyle H_{ij}=\operatorname{parity}(\texttt{bin}(i)\oplus\texttt{bin}(j))$ $$

where $\texttt{bin}(i)$ is the $m$ bit binary representation of $i$ and $\oplus$ is the bit-wise xor operation. We then define a Persian rug matrix as

$$ $\displaystyle R=n_{d}^{-1}\sum_{k\in S}H_{ik}H_{kj}$ $$

where $S\subset\{1,...,n_{s}\}$ with $|S|=n_{d}$. In fig. 1 we plot such a matrix for $n_{s}=256$, $n_{d}=40$, and $S$ chosen randomly. The matrix $R$ has diagonals equal to 1 because each diagonal is the average of $n_{d}$ terms equal to $(\pm 1)^{2}$, and a projector because the rows of $H$ are orthogonal, and thus saturates eq. 6. Furthermore, one can readily check numerically that it satisfies the statistical symmetries. There remain two variables to optimize, $a$ and $b$ (recall $\mu$ can be absorbed into $b$). We do this numerically and compare to a trained model in fig. 2.

### 4.3 Loss scaling at high sparsity

Having obtained a simple expression for the loss in terms of constants $a,b$ and two simple random variables $x\sim\text{Uniform}[0,1]$ and $\nu\sim\mathcal{N}(\mu,\sigma)$, as well has having deduced an achievable lower bound for $\sigma$, we are now able to explain why the simple ReLU model performs so well at high sparsity. For ease of notation let us use $r=\mathchoice{\frac{n_{d}}{n_{s}}}{n_{d}/n_{s}}{n_{d}/n_{s}}{n_{d}/n_{s}}$.

#### 4.3.1 Initial loss (ratio=0)

Let us first consider the $r\rightarrow 0$ limit with all other parameters fixed. Then $\sigma\rightarrow\infty$ because of the bound in eq. 6 so the fluctuations in $\nu$ overwhelms the signal term. This means that the optimal $a$ is

$$ $a=p\frac{\mathbb{E}_{u,\nu}\left[u\operatorname{ReLU}(\nu+b)\right]}{\mathbb{E }_{\nu}\left[\operatorname{ReLU}(\nu+b)^{2}\right]}+O(\sigma^{-1}).$ (8) $$

The loss then becomes

$$ $\displaystyle L$ $\displaystyle=(1-p)a^{2}\mathbb{E}_{\nu}[\operatorname{ReLU}\left(\nu+b\right) ^{2}]+p\mathbb{E}_{u,\nu}[(u-a\operatorname{ReLU}\left(\nu+b)\right)^{2}]+O( \sigma^{-1})$ $\displaystyle=a^{2}\mathbb{E}_{\nu}[\operatorname{ReLU}\left(\nu+b\right)^{2}] -2ap\mathbb{E}_{u,\nu}[u\operatorname{ReLU}(\nu+b)]+p\mathbb{E}_{u}[u^{2}]+O( \sigma^{-1})$ $$

plugging in $a$ explicitly gives

$$ $\displaystyle L=p\mathbb{E}_{u}[u^{2}]-p^{2}\frac{\left(\mathbb{E}_{u}\left[u \right]\right)^{2}\left(\mathbb{E}_{\nu}\left[\operatorname{ReLU}(\nu+b)\right ]\right)^{2}}{\mathbb{E}_{\nu}\left[\operatorname{ReLU}(\nu+b)^{2}\right]}+O( \sigma^{-1}).$ $$

Thus we can conclude that

$$ $\displaystyle\lim_{r\rightarrow 0}L=p\mathbb{E}_{u}[u^{2}]+O(p^{2})=\frac{p}{3 }+O(p^{2})$ $$

Thus we see that in the $p\ll 1$ regime we have $L\rightarrow L_{0}(p)\sim O(p)$ independent of the other parameters. We will now see that increasing $r$ will quickly cause the loss to drop to $O(p^{2})$.

#### 4.3.2 Loss upper bound

We now show that the loss drops off quickly in the sense that for $\frac{r}{p}\gg 1$ we get that $L(p)/p\rightarrow 0$, i.e. $L(p)$ scales super-linearly with $p$. We will consider the regime where $r\ll 1$ holds(^2^22For example $r=p^{1-\epsilon}$ for any $\epsilon\in(0,1)$.) so that we may take $\sigma^{2}\sim\frac{p}{r}\ll 1$.

To obtain the upper bound we will make educated estimates for values of $a$ and $\mu$ that are near optimal. In particular, in appendix A we show that the optimal value of $a$ is (after absorbing $b$ into the mean of $\nu$, $\mu$) :

$$ $a_{\text{opt}}=\frac{\mathbb{E}_{x,\nu}\left[x\operatorname{ReLU}(x+\nu)\right ]}{\mathbb{E}_{x,\nu}\left[\operatorname{ReLU}(x+\nu)^{2}\right]}.$ (9) $$

From the form of the loss, we know that $\mu$ must decrease as $p$ decreases for the loss to go down faster than $O(p)$. Thus $\nu$ has both a mean and variance approaching 0, and $a_{\text{opt}}\rightarrow 1$. Thus we plug in $a=1$ before taking these limits in the expectation of getting a good upper bound. The loss then takes the form

$$ $\displaystyle L=(1-p)L_{\text{off}}+pL_{\text{on}}$ $$

with

$$ $\displaystyle L_{\text{off}}$ $\displaystyle=\mathbb{E}\left[\operatorname{ReLU}\left(\nu\right)^{2}\right], \text{ and}$ $\displaystyle L_{\text{on}}$ $\displaystyle=\mathbb{E}\left[(u-\operatorname{ReLU}\left(u+\nu)\right)^{2}\right]$ $$

##### Off term:

The off term can be upper bounded via

$$ $\displaystyle\mathbb{E}\left[\operatorname{ReLU}(\nu)^{2}\right]$ $\displaystyle=\int_{0}^{\infty}d\nu\frac{e^{-\frac{(\nu+|\mu|)^{2}}{2\sigma^{2 }}}}{\sqrt{2\pi\sigma^{2}}}\nu^{2}=\sigma^{2}\int_{0}^{\infty}d\nu\frac{e^{- \frac{(\nu+\frac{|\mu|}{\sigma})^{2}}{2}}}{\sqrt{2\pi}}\nu^{2}$ (10) $\displaystyle\leq\frac{\sigma^{2}}{2}e^{-\frac{\mu^{2}}{2\sigma^{2}}}$ (11) $$

and thus we see we need to set $\frac{\mu}{\sigma}\gg 1$ to get a good bound. In particular, we know empirically that the loss drop happens at increasingly smaller $r$. To ensure this we let $\sigma^{2}\sim\frac{p}{r}$ scale at some rate slower than $p$. Thus to ensure that the total loss decreases faster than $O(p)$, we need $e^{-\frac{\mu^{2}}{2\sigma^{2}}}\sim O(p)$ or in other words

$$ $\displaystyle\mu\sim\sigma\sqrt{\log{\frac{1}{p}}}.$ (12) $$

##### On term:

We perform a similar, but slightly more involved computation in appendix B and combine with the off term to obtain

$$ $\displaystyle L<(1-p)O(\sigma^{2}e^{-\frac{\mu^{2}}{2\sigma^{2}}})+pO(\mu^{3}+ e^{-\frac{|\mu|^{2}}{2\sigma^{2}}}+\mu e^{-\frac{|\mu|^{2}}{2\sigma^{2}}}+ \sigma^{2}+\sigma\mu+\mu^{2}).$ $$

Plugging in the $\mu$ scaling from eq. 12 and keeping only the lower order terms gives

$$ $\displaystyle L<O\left(\sigma^{2}p\log{\frac{1}{p}}\right)\sim O\left(\frac{p^ {2}}{r}\log{\frac{1}{p}}\right).$ (13) $$

#### 4.3.3 Loss lower bound

In appendix D we also derive a lower bound in the high-sparsity limit $L>O(\frac{p^{2}}{r})$
in the high sparsity limit up to logarithmic corrections. We show this in fact holds for a more general class of activation functions. In particular, any function which acts element-wise or filters out elements will give an on-loss contribution of the form

$$ $\displaystyle\mathbb{E}_{u,\nu}[(u-f\left(u+\nu)\right)^{2}]$ $$

which has a lower bound due to $\nu$ destroying information about $u$. Thus we can conclude that

$$ $\displaystyle L\sim O\left(\frac{p^{2}}{r}\right)$ $$

up to logarithmic factors whenever $\frac{p}{r}\ll 1$.

## 5 Discussion

### 5.1 Summary

We have performed an exhaustive analysis of the toy model of superposition for permutation symmetric data in the limit of large number of sparse signals. We empirically showed that the network reliably learns statistically permutation symmetric algorithms, in the sense that the information relevant to the algorithm is encoded in the summary statistics of the elements in each row of the $W$ matrix.

Plugging these symmetries back into the reconstruction-loss allowed us to reduce its optimization problem over all model parameters to an optimization problem over three scalar variables. Of these parameters, only the noise parameter $\sigma$ was constrained by the hyper-parameters $p,n_{s},$ and $n_{d}$. We showed that $\sigma$ is minimized for a statistically symmetric $W$ that is proportional to a projector.

We then forward-engineered weights optimizing these parameters, giving the Persian rug model, and showed empirically that its loss matched the model loss closely. This shows that the trained models learn the optimal symmetric algorithm. The additional presence of small scale structure in the Persian rug model highlights the algorithm’s independence on such smaller scale details.

Finally, we considered the analytic form of the loss in the high-sparsity limit, and derived its scalings as a function of sparsity and compression ratio to be $O(p)$ when $\frac{r}{p}\ll 1$ and $\sim O(p^{2}/r)$ when $\frac{p}{r}\ll 1$, or in other words the loss drops of $O(p)$ in a ratio window scaling faster than $O(p)$.

### 5.2 Relationship to other works

#### 5.2.1 Mechanistic Interpretability and Sparse Autoencoders

Mechanistic interpretability is a research agenda which aims to understand learned model algorithms through studying their weights, (see Olah et al. () for an introduction). Recent results relating to language models include Meng et al. (), which finds a correspondence between specific facts and feature weights, along with Olsson et al. (), which shows that transformers learn in context learning through the mechanism of “induction heads”.

Elhage et al. () introduced the model we study in this paper. While that work focused on mapping empirically behaviors of the model in multiple regimes of interest such as correlated inputs, we focused on a regime with enough symmetry to solve the model analytically given observed symmetries in trained models. Chen et al. () study this model in the context of singular learning theory. As part of their work, they characterize the loss using a different high sparsity approximation than the one we present in this paper (they assume exactly one input feature is on). Then they derive a subset of the critical points and their corresponding local learning coefficients under the assumption $n_{d}=2$.

While the concept of dictionary learning was introduced by (Mallat & Zhang, ), the practical use of sparse autoencoders to understand large language models has accelerated recently due to mezzo-scale open weight models (Gao et al., ; Lieberum et al., ) and large-scale open-output models Bricken et al. (). These features are highly interpretable (Cunningham et al., ) and scale predictably. Interestingly, the scaling is quite similar for the various different architectures they consider, differing primarily by a constant, which fits with the predictions in this work.

We have seen that the dominant source of error is not from determining which features are present, but rather the actual values of those features. Small modifications to the activation functions, such as gating Rajamanoharan et al. (), k-sparse Makhzani & Frey (), or TRec non-linearity Taggart (); Konda et al. (), are insufficient to fix this problem as they do not solve the basic issue of noisy outputs. In this context our work implies that innovative architectures, that are suitable both for gradient-based training and also for decoding sparse features, must be developed.

#### 5.2.2 Compressed Sensing and Statistical Physics

It is known that compressed sparse data can be exactly reconstructed by solving a convex problem (Candes & Tao, ; Candes et al., ; Donoho & Elad, ; Donoho, ) given knowledge of the compression matrix. Furthermore, using tools from statistical physics it is possible to show that this holds for typical compressed sparse data (Ganguli & Sompolinsky, ). Learning the compression matrix is also easy in certain circumstances(Sakata & Kabashima, ). For a more general review on compressed sensing and it’s history consider the introduction by Davenport et al. (). The reconstruction procedure typically used in compressed sensing is optimizing a (convex) relaxation of finding the sparsest set of features which reproduces your data vector. This is significantly different to the setting of sparse autoencoders which try to obtain the sparse features using only one linear + activation layer.

The discrepancy between the ability of convex optimization techniques to achieve zero loss while a linear + ReLU model necessarily incurs an error suggests that a more complex model architecture is needed for sparse autoencoders when it is desirable to calculate the feature magnitude to high precision. This may occur, for example, if one wishes to insert a sparse autoencoder into a model without corrupting its downstream outputs.

#### 5.2.3 Linear Models

The linear version of the toy model we studied has garnered significant attention. Plaut () demonstrates that the minimizers of eq. 2 can be directly related to Principal Component Analysis (PCA), a fundamental technique in dimensionality reduction and data analysis (see also Bishop & Nasrabadi () for a comprehensive treatment). Kunin et al. () studies the learning dynamics of linear autoencoders when regularization is added.

### 5.3 Remaining questions and future directions

Correlated Data
We expect real world data to display significant correlations between activated features. Nevertheless, we expect there are more realistic correlation structures that remain simple enough to analyze the loss analytically. One such possibility is to use data whose correlations are described by a sparse graph.

Compressed Computation
We have shown that features can only be recovered using linear + element-wise or selection type activations together with a small error. How then are neural networks able to compute on information stored in superposition using just linear + ReLU layers? Noise injected into data early on in a computation tends to be compounded as the computation grows in circuit complexity. Thus we believe the neural networks must be at least one of the following:

- 1.
performing computation on superposed information without first localizing it,
- 2.
performing computation on superposed information by first localizing it using multiple layers,
- 3.
performing only low complexity computation on superposed information, or
- 4.
implementing a fault tolerance scheme (computation which actively corrects errors).

It would be interesting to experimentally look for these phenomena and their structure.

## Acknowledgements

We are grateful for the valuable feedback we received from Daniel Kunin, David Barmherzig, Henrik Marklund, Surya Ganguli, and Zach Furman. We acknowledge Adam Brown in particular for insightful remarks that helped us construct the Persian rug matrix. Aditya Cowsik acknowledges funding from the Simons Foundation (Award ID: 560571).

## Appendix A Optimizing over a 𝑎 a italic_a

We may optimize over $a$ analytically because it appears almost quadratically in the loss. Consider the expression for the average loss from eq. 5 with $b$ replaced with $a\cdot b$. As long as $a\neq 0$ this redefinition doesn’t change the set of accessible models.

Furtheremore let us restrict to positive $a$ which allows us to rewrite the loss as

$$ $L=\mathbb{E}_{x,\nu}[\left(x-\operatorname{ReLU}\left(a(x+\nu+b)\right))^{2}]= \mathbb{E}_{x,\nu}[(x-a\operatorname{ReLU}\left(x+\nu+b\right)\right)^{2}].$ (14) $$

The restriction to positive $a$ is acceptable because we never see negative off-diagonal elements in our trained models. Now, optimizing over $a$ is exactly linear regression; we can obtain the optimal value of $a$ with the standard method

$$ $\displaystyle 0=\frac{d}{da}L=-2\mathbb{E}\left[(x-a\operatorname{ReLU}(x+\nu+ b))\operatorname{ReLU}(x+\nu+b)\right]$ (15) $$

which implies that

$$ $a_{\text{opt}}=\frac{\mathbb{E}_{x,\nu}\left[x\operatorname{ReLU}(x+\nu+b) \right]}{\mathbb{E}_{x,\nu}\left[\operatorname{ReLU}(x+\nu+b)^{2}\right]}.$ (16) $$

Notice that the optimal $a$ is always positive, which is consistent with the assumption we made earlier.

## Appendix B Bounding the Reconstruction Error

On term:

We can start to write the on term similarly as

$$ $\displaystyle\langle(u-\operatorname{ReLU}\left(u+\nu)\right)^{2}\rangle=\left \langle\sigma^{2}\int_{-\infty}^{\infty}d\nu\frac{e^{-\frac{(\nu+\frac{|\mu|}{ \sigma})^{2}}{2}}}{\sqrt{2\pi}}(u-\operatorname{ReLU}(u+\sigma\nu))^{2}\right \rangle_{u}.$ $$

Now we write the integral in two parts to get rid of the ReLU: one when $u+\sigma\nu<0$ and one when $u+\sigma\nu>0.$ This gives

$$ $\displaystyle\underbrace{\left\langle\int_{-\infty}^{-\frac{u}{\sigma}}d\nu \frac{e^{-\frac{(\nu+\frac{|\mu|}{\sigma})^{2}}{2}}}{\sqrt{2\pi}}u^{2}\right \rangle_{u}}_{E_{r}}+\underbrace{\left\langle\sigma^{2}\int_{-\frac{u}{\sigma} }^{\infty}d\nu\frac{e^{-\frac{(\nu+\frac{|\mu|}{\sigma})^{2}}{2}}}{\sqrt{2\pi} }\nu^{2}\right\rangle_{u}}_{E_{\nu}}.$ $$

Where the first term $E_{r}$ represents error coming from the $\operatorname{ReLU}$ and the second term $E_{\nu}$ represents error coming from the noise. The scaling of $E_{\nu}$ can be easily bounded:

$$ $\displaystyle E_{\nu}<\sigma^{2}\int_{-\infty}^{\infty}d\nu\frac{e^{-\frac{( \nu+\frac{|\mu|}{\sigma})^{2}}{2}}}{\sqrt{2\pi}}\nu^{2}\sim O(\sigma^{2}+ \sigma\mu+\mu^{2})$ $$

And thus we see, unsurprisingly, that we need to set $\mu\ll 1$ to get a good bound.

To upper bound $E_{r}$ write the $u$ integral in two intervals: $[0,2|\mu|]$ and $[2|\mu|,1]$, corresponding to regions in which the interval of the $\nu$ integral does and ”decisively” does not include the the mean respectively. In particular, we have

$$ $\displaystyle E_{r}<\underbrace{\int_{0}^{2|\mu|}u^{2}du\int_{-\infty}^{-\frac {u}{\sigma}}d\nu\frac{e^{-\frac{(\nu+\frac{|\mu|}{\sigma})^{2}}{2}}}{\sqrt{2 \pi}}}_{E_{r}^{\text{mean}}}+\underbrace{\int_{2|\mu|}^{1}u^{2}du\int_{-\infty }^{-\frac{u}{\sigma}}d\nu\frac{e^{-\frac{(\nu+\frac{|\mu|}{\sigma})^{2}}{2}}}{ \sqrt{2\pi}}}_{E_{r}^{\text{tail}}}.$ $$

Since in $E_{r}^{\text{mean}}$ the $\nu$ integrals’ interval includes the mean we may as well extend the interval to the full real line to get the bound, giving

$$ $\displaystyle E_{r}^{\text{mean}}<\int_{0}^{2|\mu|}u^{2}du=O(\mu^{3}).$ $$

$E_{r}^{\text{tail}}$ can be upper bounded by setting the $\nu$ range to the maximum value of $2|\mu|$, so we have

$$ $\displaystyle E_{r}^{\text{tail}}<(1-2|\mu|)\int_{-\infty}^{-\frac{2|\mu|}{ \sigma}}d\nu\frac{e^{-\frac{(\nu+\frac{|\mu|}{\sigma})^{2}}{2}}}{\sqrt{2\pi}}= (1-2|\mu|)\int_{-\infty}^{0}d\nu\frac{e^{-\frac{(\nu-\frac{|\mu|}{\sigma})^{2} }{2}}}{\sqrt{2\pi}}<O(1)(1-2|\mu|)e^{-\frac{|\mu|^{2}}{2\sigma^{2}}}.$ $$

Putting it all together gives

$$ $\displaystyle L<(1-p)O(\sigma^{2}e^{-\frac{\mu^{2}}{2\sigma^{2}}})+pO(\mu^{3}+ e^{-\frac{|\mu|^{2}}{2\sigma^{2}}}+\mu e^{-\frac{|\mu|^{2}}{2\sigma^{2}}}+ \sigma^{2}+\sigma\mu+\mu^{2}).$ $$

Plugging in the $\mu$ scaling from eq. 12 and keeping only the lowest order terms gives

$$ $\displaystyle L<O(\sigma^{2}p\log{\frac{1}{p}})\sim O(\frac{p^{2}}{r}\log{ \frac{1}{p}}).$ (17) $$

## Appendix C Minimal Variance Bound

We will show a minimum variance bound for matrices $W$ which have all diagonals equal to 1 and also have maximum rank $n_{d}$. In this case we know that $\operatorname{Tr}W=n_{s}$. On the other hand we also know that the trace is the sum of the eigenvalues, and because $W$ has rank at most $n_{d}$ that

$$ $n_{s}=\sum_{i=1}^{n_{d}}\lambda_{i}$ (18) $$

for the eigenvalues $\lambda_{i}$ of $W$. Now we solve for the mean of the variance across rows,

$$ $\frac{1}{n_{s}}\sum_{i=1}^{n_{s}}\operatorname{Var}(\nu_{i})=\frac{4p-3p^{2}}{ 12n_{s}}\sum_{i=1}^{n_{s}}\sum_{j=1,j\neq i}^{n_{s}}W_{ij}^{2}=\frac{4p-3p^{2} }{12n_{s}}\left(\operatorname{Tr}(WW^{\dagger})-n_{s}\right).$ (19) $$

Here the first equality arises from the definition of $\nu_{i}$ (remembering that we have set the diagonals to 1 exactly) and substituting the variance of $x_{j}$, while the second equality follows because $\operatorname{Tr}(WW^{\dagger})$ is the sum of the square of all entries of $W$, and we subtract off the diagonal entries.

Because we want a bound on this quantity related to the eigenvalues of $W$, it is convenient to use the Schur decomposition of $W=QUQ^{\dagger}$. Here $Q$ is a unitary matrix and $U$ is upper-triangular with the eigenvalues of $W$ on the diagonal. This allows us to lower bound the trace

$$ $\operatorname{Tr}(WW^{\dagger})=\operatorname{Tr}(QUQ^{\dagger}QU^{\dagger}Q^{ \dagger})=\operatorname{Tr}(UU^{\dagger})=\sum_{i,j=1}^{n_{s}}|U_{ij}|^{2}\geq \sum_{i=1}^{n_{d}}|\lambda_{i}|^{2}\geq\frac{n_{s}^{2}}{n_{d}}$ (20) $$

where the last inequality follows from Cauchy-Schwarz and eq. 18. With this we find a bound on the variance

$$ $\frac{1}{n_{s}}\sum_{i=1}^{n_{s}}\operatorname{Var}(\nu_{i})\geq\frac{4p-3p^{2 }}{12}\left(\frac{n_{s}}{n_{d}}-1\right),$ (21) $$

with equality if $W$ is symmetric with all non-zero eigenvalues equal. These two conditions follow because the two inequalities in the proof become equalities when these conditions are met. This naturally leads to a candidate for the optimal choice of $W$, namely matrices of the form

$$ $W\propto OPO^{T}\text{ and }W_{ii}=1$ (22) $$

where $O$ is an orthogonal matrix and $P$ is any rank-$n_{d}$ projection matrix. This kind of matrix saturates both bounds because it is symmetric and has all nonzero eigenvalues equal to 1.

## Appendix D Lower bound on loss scaling

We now show a lower bound on the loss in the $p\rightarrow 0$ limit. To do this, we will show a more general lower bound on the on loss for any deterministic function of the pre-activation. Specifically, we would like to lower bound

$$ $\displaystyle L\leq pL_{\text{on}}$ $\displaystyle=\mathbb{E}[(u-f\left(u+\nu)\right)^{2}]$ $$

for any function $f$ with $u\sim\text{Uniform}[0,1]$ and $\nu\sim\mathcal{N}(0,\sigma)$. Recall that there is no need to consider the bias $b$ as it can be absorbed into $\nu$. Recall that the optimal function $f$ is given by

$$ $\displaystyle f^{*}(u+\nu)=\mathbb{E}[u|u+\nu].$ $$

Let’s make a change of variable from $u,\nu$ to $y\equiv u+\nu,u$, and then use the tower property to rewrite $L_{\text{on}}$ as

$$ $L_{\text{on}}=\mathbb{E}_{y\sim P_{u+\nu}}\left[\mathbb{E}_{u|y}\left[\left(u- f^{*}(u+\nu)\right)^{2}\right]\right].$ (23) $$

We first draw $y$ from the marginal distribution of $u+\nu$ and then draw $u$ from the conditional distribution given $y$. Because $f^{*}$ is exactly the conditional expectation the interior expectation becomes the conditional variance

$$ $L_{\text{on}}=\mathbb{E}_{y}\left[\operatorname{Var}\left[u|y\right]\right].$ (24) $$

Because we want to lower bound $L_{\text{on}}$ it will be convenient to start with a lower bound for the conditional variance. We will lower bound the conditional variance for $y\in[\sigma,1-\sigma]$, and then use that lower bound to find a lower bound for the loss, with a goal of showing that the loss is lower bounded by a constant multiple of $\sigma^{2}$, for $\sigma<\frac{1}{4}$. This will show that the overall loss of any strategy, even one which can perfectly estimate which features are on or off, is incapable of achieving a reconstruction error better than $O(p^{2}/r)$.

The conditional distribution for $u$ is a truncated Gaussian distribution. By Bayes’ theorem

$$ $\displaystyle P[u|u+\nu=y]$ $\displaystyle=\frac{P[u+\nu=y]P[u]}{P[y]}$ (25) $\displaystyle=\begin{cases}\frac{e^{-(u-y)^{2}/2\sigma^{2}}}{\int_{0}^{1}dxe^{ -(x-y)^{2}/2\sigma^{2}}}&\text{if }u\in[0,1]\\ 0&\text{otherwise},\end{cases}$ (26) $$

with normalizing constant $Z(y)=\int_{0}^{1}dxe^{-(x-y)^{2}/2\sigma^{2}}<\sqrt{2\pi\sigma^{2}}$. This is a truncated Gaussian distribution. Fix $y\in[\sigma,1-\sigma]$ so that all distributions are implicitly conditioned on $y$ for now. Sample $u$ via the following procedure. First we decide if $|u-y|\leq\sigma$ and then we either sample from the conditional distribution $P[u|y\text{ and }|u-y|\leq\sigma]$ or $P[u|y\text{ and }|u-y|\geq\sigma]$ with their corresponding probabilities. Let $R$ be the indicator random variable denoting $|u-y|\leq\sigma$. Then by the law of total variance

$$ $\displaystyle\operatorname{Var}\left[u~{}|~{}y\right]$ $\displaystyle=P[R=1]\operatorname{Var}\left[u~{}|~{}R=1\right]+P[R=0] \operatorname{Var}\left[u~{}|~{}R=0\right]+\operatorname{Var}_{R}\left[\mathbb {E}[u~{}|~{}R]\right]$ (27) $\displaystyle\geq P[R=1]\operatorname{Var}\left[u~{}|~{}R=1\right]$ (28) $$

where we have dropped the latter two positive terms to derive the lower bound. $P[R=1]\geq\operatorname{erf}(2^{-1/2})$ because the chance a truncated Gaussian is within one $\sigma$ of its mode is larger than that for an untruncated Gaussian, given that the truncation is more than $\sigma$ away from the mode. This condition is satisfied by construction because we have chosen $y$ to be more than $\sigma$ from the boundary.

Additionally a trivial scaling argument shows that the variance is proportional to $\sigma^{2}$ which means that there is some constant, $C>0$ such that

$$ $\operatorname{Var}\left[u~{}|~{}y\right]\geq C\sigma^{2}$ (29) $$

when $y\in[\sigma,1-\sigma]$. To complete the argument we now return to

$$ $L_{\text{on}}=\mathbb{E}_{y}\left[\operatorname{Var}\left[u|y\right]\right] \geq\mathbb{E}_{y}\left[\operatorname{Var}\left[u|y\right]1_{y\in[\sigma,1- \sigma]}\right]\geq C\sigma^{2}P[y\in[\sigma,1-\sigma]].$ (30) $$

For $\sigma=1/4$ this probability is clearly finite and for $\sigma<1/4$ it is increasing as $\sigma$ decreases so it is uniformly bounded below by a constant $C_{1}$. So finally

$$ $L_{\text{on}}\geq C^{\prime}\sigma^{2}\implies L\geq C^{\prime}p\sigma^{2} \approx\frac{C^{\prime}p^{2}}{r}.$ (31) $$

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. How does the solution proposed by the authors differ from the dictionary learning solution? Can the authors delineate conditions under which trained autoencoders converge to their solution instead of overcomplete, incoherent dictionaries?

2. The previous literature does not study autoencoders and sparse coding in the context of interpretability. Are there specific questions that arise in this context that require the authors' model? Does the authors' analysis answer these questions?
soundness: 2
strengths: The empirical observations and theoretical results are reasonable and the authors provide an adequate description of the solution $W$.
summary: This paper characterizes the solutions of training sparse autoencoders (defined as $f(\mathbf{x}) = \textrm{ReLU}(W_{out} W_{in} \mathbf{x} + \mathbf{b})$) with the reconstruction loss $L=\mathbb{E} || \mathbf{x} - f(\mathbf{x}) ||^2$. The paper first trains an autoencoder on synthetic data consisting of sparse vectors generated from a Bernoulli-Uniform model and records some empirical observations about the composite matrix $W=W_{out} W_{in}$. The authors note that the off diagonal entries of this matrix seem to follow a normal distribution that seems to be identical across rows and that the diagonal entries are also identical. This allows them to reduce the loss function to a function of three scalar variables corresponding to the diagonal, off-diagonal, and bias terms. The solution they obtain is of the form $W = OPO^\top$ where $O$ is orthogonal and $P$ is a rank $n_d$ projection matrix. The authors also characterize how the loss varies as a function of the Bernoulli sparsity parameter $p$.
weaknesses: Connections between sparse coding and autoencoders have been explored extensively in the literature. This paper does not cite or engage with any of this literature. 

The empirical observations about $W$ having identical diagonal entries and the off diagonal entries being small and having a common distribution can be captured by $W = A^\top A$ where $A \in \mathbb{R}^{n_d \times n_s}$ is an overcomplete, incoherent dictionary with unit norm columns. This ensures that the diagonal terms are identical while the off-diagonal terms are small in magnitude. Moreover it is well known that random matrices satisfy the incoherence/restricted isometry condition and a random $A$ ensures that the off-diagonal entries follow the same distribution.

This model of dictionary learning, along with Bernoulli-Uniform or Bernoulli-Gaussian sparse codes has been studied in the complete and overcomplete setting with neural and non-neural algorithms [3,4,5].

Reference [1] studies how ReLU autoencoders with tied weights can learn dictionaries with data from this generative model. The authors show (in theorem 3.1) that one layer of a ReLU autoencoder can recover the support of the sparse code with high probability if the weights are close to the ground truth dictionary. The result also proposes a bias vector with all negative entries whose magnitude can be derived from the incoherence of the dictionary. The paper also uses a landscape argument to show that the ground truth dictionary is a critical point for the reconstruction loss. Reference [2] goes further to show that autoencoders trained to minimize reconstruction error with gradient descent can recover the ground truth dictionary.

[1] Rangamani, A., Mukherjee, A., Basu, A., Arora, A., Ganapathi, T., Chin, S., & Tran, T. D. (2018, June). Sparse coding and autoencoders. In 2018 IEEE International Symposium on Information Theory (ISIT) (pp. 36-40). IEEE.

[2] Nguyen, T. V., Wong, R. K., & Hegde, C. (2019, April). On the dynamics of gradient descent for autoencoders. In The 22nd International Conference on Artificial Intelligence and Statistics (pp. 2858-2867). PMLR.

[3] Arora, S., Ge, R., Ma, T., & Moitra, A. (2015, June). Simple, efficient, and neural algorithms for sparse coding. In Conference on learning theory (pp. 113-149). PMLR.

[4] Agarwal, A., Anandkumar, A., Jain, P., & Netrapalli, P. (2016). Learning sparsely used overcomplete dictionaries via alternating minimization. SIAM Journal on Optimization, 26(4), 2775-2799.

[5] Spielman, D. A., Wang, H., & Wright, J. (2012, June). Exact recovery of sparsely-used dictionaries. In Conference on Learning Theory (pp. 37-1). JMLR Workshop and Conference Proceedings.

[6] Refinetti, M., & Goldt, S. (2022, June). The dynamics of representation learning in shallow, non-linear autoencoders. In International Conference on Machine Learning (pp. 18499-18519). PMLR.

While these prior results study a slightly different model architecture, they can be adapted to the model studied in this paper with a little effort. In my opinion the results obtained by the current paper can be explained by the sparse coding model, and it is not true that this has not been studied in the context of autoencoders.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: My questions are mostly overlapped with the weakness section. I appreciate it if the authors may further instantiate the following points:

W1. What is the interpretation of the model if the input is not sparse? Is the symmetry property and the relevant conclusion partially hold? Does this mean the interpretation of the weights and the optimal W will change significantly if the compression ration becomes high?  

W2. Generally speaking, what is the significance of the proposed model if the input is no longer sparse, whereas in most cases the input is not sparse?
soundness: 3
strengths: S1. Novelty: The paper is novel in it both theoretically and empirically performed exhaustive analysis of autoencoder trained on sparse data, interpreting the essential structure needed to represent such data. I like the way they demonstrates the distribution of the weights. 

S2. The paper firstly reveals that a statistically symmetric strategy exists in certain regimes of the macroscopic parameters and thus optimal choice for W exisits based on the proposed  Persian Rug Model. The theoretical support of the motivation is strong. 

S3. The presentation of the paper is good, with a good clarity.
summary: The paper both theoretically and empirically performed an exhaustive analysis of autoencoder trained on sparse data,
extracting the essential structure needed to represent such data. The paper firstly reveals that the change of loss as the compression ratio varies and shows that when compression ratio is near zero, the loss drops to zero entirely. This motivates the authors to consider the structure of the weight in both linear and non-linear cases. They artificially include the statistical permutation symmetry into the weight and found the artificial weight and the trained model returned the same result. The conclusion is they qualitatively show that a statistically symmetric strategy exists in certain regimes of the macroscopic parameters. Thus optimal choice for W exisits based on the proposed  Persian Rug Model.
weaknesses: W1. What is the interpretation of the model if the input is not sparse? Is the symmetry property and the relevant conclusion partially hold? I noticed the theories are mostly based on the compression ratio. Does this mean the interpretation will change if the compression ration becomes high?  

W2. What is the significance of the model if the input is no longer sparse?

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 4
questions: - (Minor) Do the authors believe the study could be extended to non-isotropic cases, where some features are of higher importance, e.g. in the block case where blocks of different features have the same importance ?

- Is it correct that (7) also corresponds to the weights a linear model would learn ? If that is the case, does the improvement in reconstruction loss from the reLU model come from the bias and activation alone ? 

- As pointed above, is there any specific reason the particular realization of 4.2.2 of optimal weights is of specific interest, compared to randomly drawing an orthogonal matrix $O$ and constructing a weight matrix as in (7)? Do the authors claim the model learns weights exactly corresponding to 4.2.2, or simply that the loss coincide ?
soundness: 4
strengths: On a technical level, the authors provide a detailed study of the optimal population loss achieved by the considered model, which is non trivial given the ReLU non-linearity. To the best of my awareness, this derivation is new compared to the original work of Elhage et al. (2022). The finding of the weight matrix achieving the same loss as the model optimized using classical gradient-based optimizers is an interesting construction. Although I have not checked the math in detail, the derivation looks scientifically sound.

On a writing level, the paper is very clearly written, with sufficient discussion and intuition building. The assumptions are illustrated by careful numerical experiments.  Sufficient discussion accompanies all important equations.
summary: The authors study the compression of sparse data using a two-layer autoencoder with ReLU output non-linearity, considering the toy model of Elhage et al (2022). Building on an ansatz for the product of the weights, supported by numerical observation, they derive a closed-form expression for the population loss depending on only three scalar parameters, corresponding to the bias, common diagonal element of the product matrix, and variance of off-diagonal elements. They deduce from it a lower-bound on the loss, and construct a deterministic weight matrix that achieves it. They find that the corresponding loss coincides with that of the trained model.
weaknesses: My main concern lies in the scope of the contribution. Some of my concerns are tied to my questions, which I list in the next section.

- To my understanding, the study is confined to the isotropic case of Elhage et al. (2022), as it builds upon the ansatz that all diagonal elements of $W_{out}W_{in}$ are equal. This means there is not really any feature learning in the model, although the considered setting remains non-trivial.

- While the study is interesting, it does not add much more phenomenological insights compared to the original work. In particular, the authors construct a matrix which achieves the same optimal population loss as the trained model, which is one of the strengths of the work. However, I understand any matrix of the form (7) is similarly optimal, and I do not understand why the explicit realization in 4.2.2 brings any further insight. I might have misunderstood a point, and am happy if the authors bring clarifications on this point.

Overall, I think the study is technically sound, the insights tied to the analysis of interest, but brings little novel phenomenological insight.  I am thus giving an accept score -- because the paper looks sound --, but not a high one.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 2
questions: 1. On line 111, why must we construct $x_i$ like this? Can we sample $u_i$ from the other distributions, like uniform[-1, 0]? This still makes $x$ sparse, right? I guess sampling $u_i$ from the other distributions may not work because you AE has ReLU at the output activation? Are you assuming that $x$ are outputs of ReLU or any other positive activation? 
    * How well does this synthetic distribution of sparse data agree with the distribution of sparse data in reality? Here you are assuming its dimensions to be i.i.d, in reality this is probably not the case.
2. You said “We must restrict to W with rank no more than nd” on Line 324. Do we need to take this rank condition into account when deriving Equation (4) and making the following reasonings from Line 332 to 348? If not, why? Could $\nu_i$ correlate with each other due to the low rank of $W$? Are you making any i.i.d. assumption about the off-diagonal elements here?
    * Also, this “because νi becomes Gaussian” seems to be a conjecture based on your experiment in Section 3.2. I’m not familiar with your test so I don’t know how reliable it is. That said, you’d better make this assumption stand out with something like `\newtheorem{asp}{Assumption}`. 
3. Given your Persian rug construction of $W$, any way to retrieve $W_\text{in}$ and $W_\text{out}$ from it? Are they unique? Can we use the pair on Line 398? 
4. Expand your derivation of Eq. 19 in Appendix B to make it clearer. Why is there no $a$ in it (considering your definition of $\nu$ on Line 333)? Does this equation only hold when you construct $x_i$ as described on Line 111?
soundness: 3
strengths: 1. Has theoretical analysis that looks legit to me, even though it is only focused on a very specific type of autoencoders on a specific type of datasets. 
2. The observation leading to theoretical analysis structure of this work can effectively convey the authors’ motivation.
summary: The work tries to derive analytical ways to obtain the near-optimal solutions of a specific type of autoencoders for synthetic sparse datasets. After theoretical analysis, it provides the process for constructing the autoencoder’s weight matrix without numerical optimization.
weaknesses: The main weakness of this paper is its inconsistency and the absence of experiments for supporting their theoretical findings in Section 4, which greatly undermines the soundness of their results. Also, the authors should avoid introducing jargons and verbosity, be concise and accurate about what they mean, or provide necessary references if the jargons are widely accepted elsewhere. 
I suggest that the authors spend more time streamlining their work to improve its readability and precision, and bolster their theoretical claims with empirical results.

Regarding contribution, the results in this paper might be too restricted to a specific type of problem and model. I'm not sure how crucial this problem is to LLM development, and because I also value rigorous theoretical studies of problems no matter how small they seem, I refrain from saying harsh words like "it's trivial". That said, the authors had better show their problem's importance with better literature review and critical discussion. 

## Major

### Abstract

1. Some claims in the abstract are not (well) reflected in the main text. For instance, “changes to the elementwise activation function or the addition of gating can at best improve its performance by a constant factor”. I guess this is about the discussion on Line 486 to 492. Experiments need to be done to verify this.

### Section 1

1. Line 76: Please provide clear definition of “permutation symmetric data” and “thermodynamic limit”. At least they do not appear in the work (Elhage 2022) you cited, and they sound kinda vague and unintuitive. For instance, how can “permutation symmetric” be derived from “no input feature is privileged”? As long as the features are not identical you cannot say it’s “permutation symmetric or invariant”, right? If these are jargons commonly used in your out-of-ML research area, you need to translate it for your mortal audiences in the ML community. 

### Section 2

1. How many samples are there in your dataset on Line 111? What is the rank of this data matrix you sampled? Could the rank of this data matrix affect or even undermine your later observations and motivations? Did you scale up the number of samples as you increase $n_s$? If so, how? You need to provide more details about your experiments.

### Section 3

1. On line 142 there is *“due to an immediate drop in the loss starting around nd/ns ≈p. The slope and duration of this initial fall is controlled by p. In particular, in the high-sparsity regime (pclose to zero), the loss drops to zero entirely near the nd/ns ≈0 regime.”* 
I don’t think any part of Fig. 2b can clearly show this, especially considering that the upper bound of the x-axis (0.8) is way larger than p=0.05. You need to plot the loss curve at different p to demonstrate that, and zoom in on the “nd/ns ≈0 regime” and “nd/ns ≈p regime” to illustrate your points. Also, use the `subfigure` environments to enclose the subfigures.
2. On line 152, there is *“nd/ns →∞ limit”*, but Fig. 3,4,5,6 cannot show any of such limit. All nd/ns are under 0.8.  Do you mean “nd/ns →0 limit” (although this makes no sense), or just ns →∞? Because of this, I'm uncertain when reading your statements from Line 155 to 161. 
    * If it's the latter, as I understand it, 8192 might still be kinda small. 
    * Are you always assuming that nd < ns throughout the entire paper? Be clear.
3. Should provide a clear definition of “permutation invariant” mathematically. Also, what is “permutation symmetric” matrix? Be clear about what you say, or provide citations at least. 
4. Elaborate on “We hypothesize this occurs because for small models (and small feature probabilities) the contribution to any output from off-diagonal terms may fluctuate because of the small size of the matrix.” on Line 215. It’s like saying many things without saying anything. 
5. Elaborate on Line 236 to 239. Actually, this sentence is too convoluted to read. 
6. On Line 236 “Similar features are present in the two plots”. On Line 240 “we give three plots”. What and where are these five plots? Be specific. 
7. For Line 266 to 269 and Fig. 5. What is $\nu$? What is $\Delta \text{var}(\nu)$? What is “the fluctuation of the total norm of the off-diagonal elements in each row in W across rows”? Is it simply the variance of the off-diagonal elements’ variances? What do you mean by signal-to-noise ratio? What is considered as noise and what is not? The writing surely needs to be improved. If the concepts are introduced in the later sections (like in (4)), you’d better considering rearranging them to improve your narrative.
8. Elaborate on your choice of $\Lambda$ in (3). What is your motivation? What does it indicate? Any related literature where this metric is used and explained in detail? What is its relation to the Lyapunov condition? You say “strong numerical evidence”, how strong? Provide introductions or references for your readers to understand. 
9. Rephrase Line 304. Unclear. Better give a rigorous math expression for clarity. 

### Section 4

1. The important variable σ is not introduced very clearly. The only related explanation is *“σ characterizing the root mean square of the off diagonal rows”*. I can hardly get what you mean. Please use a better expression. 
2. You need to show your Persian rug’s optimality with experiments. 
3. Does Eq. 8 imply we can fix b at any value and only adjust a to obtain the optimal algorithm? If not, how to select $b$?
4. Verify "This implies that even if ... scaling of the loss for $r>p$." on Line 490 with experiments.

## Minor

1. Always cite related works when you introduce new concepts, like Hadamard matrix, Lyapunov condition, etc.
2. Unnecessary verbosity. For instance:
    1. (Line 185 to 188) “…diagonals are large…as possible.” Why don’t you just say $W \approx c I$?
    2. (Line 22 to 25) “Unlike…oblivious”. Why should the algorithm care?
    
    You should flesh out your paper with rigorous and detailed definitions of your important concepts rather than these.
    
3. What do you mean by “r=0.25” in Fig. 2a? Is r defined on Line 424? 
4. Line 178 “p = 4.5% and ratio ns = 512.” Is it a typo? Do you mean nd/ns = 512?
5. On line 204, please provide a clear definition of “mean-square fluctuation” or $\Delta$. Is it simply the variance (looks so on Line 212)? In addition, variance is not sufficient for evaluating if all elements in a set are nearly the same. You’d better additionally evaluate the maximum absolute difference. I also suggest replacing variance with standard deviation.
6. Why use $A$ rather than $a$ in Appendix A? Improve consistency.
8. What is $\mathbf{B}$ on Line 323 and 363? $\mathbf{b}$ in (1) and (2)? Improve consistency.
9. Use the correct citation command: \citep and \citet.
10. Use “Fig. xx” instead of “fig. xx” to refer to figures, and "Eq. xx" instead of "eq. xx" for equations.

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: See the weaknesses section.
soundness: 2
strengths: The paper provides to some extent a solid analysis of the proposed model.
summary: The paper studies the following auto-encoder model: linear encoding + linear decoding with additional denoising (shifted-ReLU) application on top. This model is applied in context of compression of i.i.d. sparse uniform data. The authors provide an analysis of the model behaviour (MSE bound) and identify the set of matrices which matches the performance of the trained model.
weaknesses: The whole interpretability premise is evidently far fetched especially whenever LLMs appear in the related context. The paper does not provide any evidence that the phenomenology observed in their i.i.d. + small model setting appears on any even "toy"-ish real data, or justify enough whether the prescribed methodology could be put forward to gain an understanding of the behaviour of the real world models. The suggestion is to reframe and focus on the subject of substance - analysis of the autoencoder model.

I fail to see why so much of the paper space is occupied by a the "Gaussian-ity" experiments in Section 3.2. The authors do no provide any rigorous training results anyway (that the model will converge to a specific subset of matrices). In this view, since the proposed behaviour is assumed anyway using the empirical evidence, the amount of histograms and the accompanying description is clearly an overshot (which does not really add any further weight to the submission).

The technical machinery used for the analysis is not of a separate interest. 

The main concern, however, is the emergence of a certain structure in the weights. In particular, it is completely unclear whether the "Persian carpet" structure is indeed favoured by the network implicitly or just a set of weights which achieves the close to optimal MSE.
If later, the observations in the paper does not shed any light on any sort of "interpretability", since the network does not implement the proposed weight pattern when trained. Prior, however, is not supported by any evidence (no SGD/gradient flow convergence result).

The phenomena of uniform enough weights is not surprising on its own given that the signal is i.i.d., especially in the view of related work.
Specifically, could authors provide any intuition why specifically "Persian carpet" (except visual appeal)? There is no evidence whatsoever that a simpler uniform design will not achieve desired. Why not use some other design: for instance, rotationally invariant matricies with specific spectra (if the input distribution is centered, see next paragraph) and exploit Approximate Message Passing analysis (granted, the optimal denoising will not be close to ReLU)?

Another concern: why the data is not centered (which is *common practice*)? It is hard to judge, in this case, how much of an effect on the proposed structure it has, in the view of ReLU - which has cut-off at zero. 

I also assume that in equation (2) there is normalization in dimension missing. Otherwise, dropping "i" in (5) does not make sense.

Lastly, the particular choice of ReLU model (1) in the view of above remains quite unclear. It is hard to judge if any of the observations are not an artifact of the modeling.

The authors fail to acknowledge *a lot* of the related literature which is non-trivially connected to the observed model behaviour. Below, a non-exhaustive list:
-  AEs: (Refinetti & Goldt, 2022); (Cui & Zdeborova, 2023); "Analysis of feature learning in weight-tied autoencoders via the mean field lens" Nguyen (2021); (Shevchenko et al., 2023); (Kögler et al., 2024);
- Approximate Message Passing algorithms: (Donoho et al., 2009); rotationally invariant design (Rangan et al., 2019; Schniter et al., 2016; Ma & Ping, 2017; Takeuchi, 2019); optimal spectral design (Ma et al., 2021).

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 2
questions: - In Figure 1, you show that the structure of the artificial weights resembles a Persian rug. Out of curiosity, what do the trained weights look like?

- On line 155, you mention that in the $r \to \infty$ limit, the diagonal elements of $\mathbf{W}$ become the same, while the off-diagonal elements become zero-mean noise, resulting in $\mathbf{W} = a \mathbf{I}$. However, since $\mathbf{W}$ is the product of two low-rank matrices, shouldn’t it have a rank no greater than $n_d$? How, then, does it approach a full-rank matrix?

- On line 355, the symbol $\mu$ is not defined. Could you also clarify the main message of this paragraph?

- On line 406, the parity function is not defined.

- On line 410, should the LHS be $R_{ij}$? It would also be helpful to include a theorem stating that $\mathbf{R}$ satisfies the statistical symmetries, with a proof deferred to the appendix.
soundness: 3
strengths: - This paper presents a novel study on the effects of input symmetry in superposition models.  
- This paper identifies a statistical permutation symmetry in the learned model weights, which is verified by designing a "Persian rug model" that achieves similar performance to the trained model.  
- The authors provide theoretical upper bounds for the reconstruction loss, which their experiments confirm by showing a notable drop in loss around $r\approx p$.  
- Detailed explanations throughout the paper improve the understanding and interpretation of the experimental results.
summary: This paper examines a toy model of superposition, revealing that when the input data exhibits permutation symmetry, the auto-encoder learns an algorithm responsive only to macroscopic parameters. Using this insight, the authors derive upper bounds on the reconstruction loss and design hand-crafted symmetric models that match the loss of the trained model.
weaknesses: - **Writing**: The empirical observations in the paper feel overly repetitive, while the theoretical sections come across as overly technical and challenging to follow. For example, Section 3 spends three pages describing that the diagonal elements and bias become uniform while off-diagonal elements become noisy. Section 4 includes excessive technical detail in the derivations, which could be condensed into formal theorems, with the technical details moved to the appendix for easier readability.

- **Theory**: Many theoretical discussions in this paper lack rigor. For instance, the discussion between lines 353 and 361 is confusing and does not seem to follow rigorous mathematical reasoning. I will provide further points of confusion in the questions section. Notably, as a theoretically-focused paper, it lacks any formal theorems, making it difficult to digest. Summarizing results into well-defined theorems with clearly explained notations would improve clarity.

- **Experiments**: The experimental scope appears limited. The paper only presents results on a toy model with input data following a sparse, uniform, permutation-symmetric distribution. However, real-world data, such as natural language, does not exhibit permutation symmetry. Including experiments on real-world data would help demonstrate the applicability of these results.

In summary, the experimental scope of this paper seems narrow, and the theoretical sections lack the rigor expected of a theoretical work. I recommend that the authors refine their presentation of theoretical results and expand their experiments. I’m open to any clarifications if there has been a misunderstanding of the paper's theory from my end.

Meta-review: The paper studies the compression of sparse data via a two-layer autoencoder with ReLU output non-linearity. Building on an ansatz for the product of the weights (supported by numerical evidence), the authors derive a closed-form expression for the population loss depending on three scalar parameters. This leads to a lower-bound on the loss, which is achieved by a certain deterministic weight matrix.

Reviewers ru6M and dnUr have raised important issues concerning the absence of comparisons with relevant related work (providing a rather extensive list of references) and the technical contributions themselves (no analysis of global optimality; little justification about the parameterization considered in the paper; unclear whether the proposed solution actually occurs in practice). Such concerns were unfortunately not resolved during the discussion phase and therefore I recommend a rejection at this stage. 

I do encourage the authors to address the comments of the reviewers in a thorough revision and submit an improved version to a future venue.

Decision: Reject