Motivation: This paper addresses representational collapse in Transformers, a fundamental issue affecting model expressiveness. While Dong et al. (2021) showed pure attention suffers rank collapse, their analysis excluded components present in real Transformers like BERT and GPT.

Mathematical setup: Consider a Transformer layer with input

\[X^{(l)} \in \mathbb{R}^{n \times d_{\text{model}}}\]

where $n$ is sequence length. Multi-head attention computes:

\[\text{MHA}(X) = \sum_{h=1}^H Y^{(h)} W_O^{(h)}\]

where

\[Y^{(h)} = A^{(h)} X W_V^{(h)} \in \mathbb{R}^{n \times d_k}\]

and

\[A^{(h)} = \text{softmax}(XW_Q^{(h)} W_K^{(h)T} X^T / \sqrt{d_k})\]

Layer normalization applies:

\[\text{LN}(x_i) = \frac{x_i - \bar{x}_i}{\sigma_i} \odot \gamma + \beta\]

The residual update is:

\[X^{(l+1)} = X^{(l)} + \text{MHA}(X^{(l)}) + \text{FFN}(\cdot)\]

Assumptions:

1. Non-degeneracy: no constant rows, $\sigma_i > 0$
2. $\mathbf{1} \notin \text{rowspace}(X)$
3. Scale parameters $\gamma_j \neq 0$

4. Standard deviations not in column space of centered data

   Toy example: When $d_{\text{model}} = 4$, $H = 2$, $d_k = 2$, if heads specialize to 2D subspaces that align, then $\text{rank}(\text{MHA}(X)) \leq 2 \ll 4$, creating a bottleneck the MLP must expand.

   Formal objective: Characterize when

   \[\text{rank}(X^{(l+1)}) > \text{rank}(X^{(l)})\]

   and quantify the head-channel non-identifiability dimension

   \[\dim(\mathcal{K}_{h^*}) = n(H-1)d_k\]

The paper provides three main theoretical results:

1. Layer Normalization Rank-Neutrality: Under non-degeneracy conditions, layer normalization preserves affine rank exactly:

\[\text{arank}(\text{LN}(X)) = \text{arank}(X)\]

2. Residual Rank Obstruction: For the update $X^{(l+1)} = X^{(l)} + \text{MHA}(X^{(l)})$, generically:

\[\text{rank}(X^{(l+1)}) > \text{rank}(X^{(l)})\]

This uses measure-theoretic arguments showing rank-preserving inputs form a measure-zero set.

3. Head-Channel Non-Identifiability: Given only the sum $\sum_h Y^{(h)}W_O^{(h)}$, individual head contributions cannot be recovered. The recovery ambiguity dimension is:

\[\dim(\mathcal{K}_{h^*}) = n(H-1)d_k\]

Applied to toy example: With $H=2$, $d_k=2$, $n=4$: the ambiguity dimension is $4 \times 1 \times 2 = 8$ degrees of freedom per layer. Even knowing the mixed output perfectly, recovering individual head contributions requires 8 additional parameters.

Prior work:

1. Dong et al. (2021): showed pure attention without residuals suffers doubly-exponential rank collapse; proposed MLP as remedy
2. Wu et al. (2024): analyzed LayerNorm’s complex, regime-dependent effects on collapse rates

3. Nait Saada et al. (2024): identified rank collapse in width (distinct from depth collapse)

   Delta: This paper adds three corrections to Dong’s framework:

   • LayerNorm is precisely rank-neutral (preserves affine rank exactly), not “plays no role”
   • Residual connections generically prevent collapse without MLP contribution
   • Identifies head-channel non-identifiability as distinct from rank collapse

   Theory-specific assessment:

   • Main theorems are somewhat predictable extensions of known linear algebra, but the precise characterization is valuable
   • Proofs use standard measure theory and determinantal varieties - routine but carefully executed
   • No comparison with lower bounds provided; tightness unclear
   • The head-channel non-identifiability phenomenon appears genuinely new

   Verdict: INCREMENTAL — solid theoretical refinements of existing collapse analysis, but the core insights about residuals preventing collapse and LayerNorm being rank-neutral follow fairly directly from standard linear algebra.

1. Rank-Neutrality Theorem: Uses step-by-step analysis of LayerNorm operations:

• Mean subtraction:

\[\tilde{X} = XP\]

where

\[P = I - d^{-1}\mathbf{1}\mathbf{1}^T\]

• Since $\mathbf{1} \notin \text{rowspace}(X)$, projection $P$ is injective on rowspace
• Row-wise rescaling by $D = \text{diag}(1/\sigma_1, \ldots, 1/\sigma_n)$ preserves rank generically
• Affine rescaling $y_i = \hat{x}_i \odot \gamma + \beta$ preserves affine rank when $\gamma_j \neq 0$

2. Residual Rank Obstruction: Uses algebraic geometry and measure theory:

\[\{B : \text{rank}(A + B) \leq \text{rank}(A)\} = V_r - A\]

where $V_r$ is the determinantal variety of matrices with rank $\leq r$. Key inequality:

\[\text{codim}(V_r) = (n-r)(d_{\text{model}}-r) \geq 1\]

Since $V_r$ has measure zero, its translate also has measure zero.

3. Head-Channel Recovery: Constructs the recovery subspace:

\[\mathcal{K}_{h^*} = \ker\left((Y^{(1)}, \ldots, Y^{(H)}) \mapsto Y^{(h^*)}W_O^{(h^*)} \mid \sum_h Y^{(h)}W_O^{(h)} = 0\right)\]

Uses dimension counting in the joint parameter space $\mathbb{R}^{nHd_k}$ mapping to $\mathbb{R}^{nd_{\text{model}}}$.

Purely theoretical with some architectural analysis of BERT-base parameters. The paper provides concrete numbers: for BERT-base with $n=512$, $H=12$, $d_k=64$, the head-channel non-identifiability dimension is 360,448 per sequence (704 per token) per layer.

The paper proposes Position-Gated Output Projection (PG-OP) as a constructive remedy with <1.6% parameter overhead, but provides no empirical validation. Empirical validation would require:

1. measuring actual rank evolution in trained Transformers
2. correlating head specialization with MLP width requirements
3. testing whether PG-OP improves head interpretability without performance loss.

Limitations:

1. Non-degeneracy assumptions H1-H4 are NATURAL (widely satisfied in practice and maintained during training)

2. Generic statements via measure theory are TECHNICAL (hold almost everywhere but give no finite-sample guarantees)

3. Analysis focuses on linear algebraic properties, not optimization dynamics - TECHNICAL limitation

4. Head-channel non-identifiability remedy (PG-OP) is only partial and adds parameters - RESTRICTIVE for deployment

5. No connection to actual task performance or model quality - TECHNICAL gap between theory and practice

   Open problems:

6. Finite-sample bounds: Convert generic (measure-theoretic) statements into finite-sample probability bounds for practical architectures

7. Complete head-channel remedy: Develop architecture modifications that fully resolve head-channel non-identifiability without significant parameter overhead or performance loss

Motivation: Continuous causal discovery methods couple representation learning with non-convex acyclicity optimization, leading to local optima and poor scalability in high dimensions. Extracting directed acyclic graphs (DAGs) from observational data remains computationally bottlenecked by optimization rather than statistical estimation.

Mathematical Setup: Consider observational data generated by a continuous Additive Noise Model (ANM). Let $G = (V, E)$ be a DAG over $d$ variables with structural equations:

\[x_i = f_i(x_{pa(i)}) + \epsilon_i, \quad i = 1, \ldots, d\]

where $pa(i)$ denotes parents of node $i$, $f_i$ are smooth functions, and $\epsilon_i$ are mutually independent noise variables.

Assumptions:

1. Homoscedastic Gaussian noise: $\epsilon_i \sim \mathcal{N}(0, \sigma^2)$ with variance $\sigma^2 > 0$
2. Causal faithfulness: For any edge $j \to i \in E$, we have $\mathbb{E}_{p(x)}[(\partial f_i/\partial x_j)^2] > 0$

3. Score function $s(x) = \nabla_x \log p(x)$ is twice differentiable

   The Score-Jacobian Information Matrix (SJIM) is defined as:

   \[I = \mathbb{E}_{p(x)}[-\nabla^2_x \log p(x)]\]

   Toy Example: When $d=2$ with linear functions $x_1 = \epsilon_1$ and $x_2 = \beta x_1 + \epsilon_2$, the SJIM becomes:

   \[I = \frac{1}{\sigma^2} \begin{pmatrix} 1 + \beta^2 & -\beta \\ -\beta & 1 \end{pmatrix}\]

   The leaf node $x_2$ has diagonal entry $I_{22} = 1/\sigma^2$, while the root has $I_{11} = (1 + \beta^2)/\sigma^2 > 1/\sigma^2$.

   Formal Objective: Extract the topological ordering $\pi$ by iteratively identifying leaf nodes via:

   \[\ell = \arg\min_i I_{ii}\]

SSTS Algorithm: The Score-Schur Topological Sort operates in two decoupled stages:

Stage 1 - Generative Modeling: Train an unconstrained score network $s_\theta(x)$ by minimizing:

\[\min_\theta \mathbb{E}[\|s_\theta - \nabla \log p\|_2^2]\]

Stage 2 - Algebraic Extraction:

1. Construct empirical SJIM: $\hat{I} = \frac{1}{2}(\hat{J} + \hat{J}^T)$ where $\hat{J} = \frac{1}{

   D

   }\sum_{x \in D} -\nabla_x s_\theta(x)$

2. Apply Group Lasso sparsity prior for high dimensions:

   \[L_{sparse} = \lambda_{sparse} \sum_{j=1}^d \|[W_1]_{:,j}\|_2\]

3. Block-SSTS marginalization: Group parallel nodes using tolerance $\gamma$:

   \[B = \{k \in S | \hat{I}_{kk} \leq \min_{i \in S}(\hat{I}_{ii}) + \gamma|\min_{i \in S}(\hat{I}_{ii})|\}\]

4. Compute Tikhonov-regularized Block Schur complement:

   \[\hat{I}_{S',S'} \leftarrow \hat{I}_{S',S'} - \hat{I}_{S',B}(\hat{I}_{B,B} + \lambda_{ridge}I)^{-1}\hat{I}_{B,S'}\]

   Application to Toy Example: For the 2D linear case, SSTS identifies node 2 as the leaf (minimum diagonal), removes it via Schur complement, leaving the 1×1 submatrix for node 1. The extracted order is $\pi = [2, 1]$, correctly capturing the causal direction $x_1 \to x_2$.

Prior Work:

1. NOTEARS (Zheng et al., 2018): Introduced continuous algebraic acyclicity constraint $\text{tr}(e^{W \circ W}) = d$ enabling gradient-based optimization
2. Score-based causal discovery (Rolland et al., 2022): Used score function Jacobians to identify leaf nodes but required recursive model retraining

3. DAGMA (Bello et al., 2022): Advanced continuous optimization with log-determinant acyclicity characterization

   Delta: This paper introduces the mathematical equivalence between iterative leaf marginalization and Schur complement operations on the SJIM. The key insight is that graph marginalization can be performed algebraically without retraining score models.

   Theory-specific Assessment:

   • Main theorem: The exact equivalence (Theorem 3.2) between marginalization and Schur complement in linear ANMs is surprising and non-obvious. The extension characterizing the expectation gap (Proposition 3.4) in non-linear cases provides new theoretical insight.
   • Proof technique: The proof leverages properties of Gaussian graphical models and Fisher information matrices in novel ways. The characterization of the expectation gap via gradient covariance is genuinely new.
   • Tightness: No lower bounds are established. The approximation quality depends on the magnitude of $\text{Cov}(\nabla f_l)$ which is problem-dependent.

   The theoretical contribution is substantial - proving exact algebraic equivalence in linear cases and precisely characterizing the approximation gap in non-linear cases represents a clear advance over prior heuristic approaches.

   Verdict: SIGNIFICANT — The mathematical mapping of causal marginalization to matrix algebra provides both theoretical insight and computational advantages that most causal discovery researchers should understand.

The proofs employ several key technical components:

1. Leaf Identifiability (Theorem 3.1): Uses the Markov property and chain rule to decompose the joint log-likelihood. The key insight is isolating the contribution of children nodes:

   \[I_{ii} = \frac{1}{\sigma^2} + \frac{1}{\sigma^2}\sum_{j \in \text{ch}(i)} \mathbb{E}\left[\left(\frac{\partial f_j}{\partial x_i}\right)^2\right]\]

   The proof exploits the independence of $\epsilon_j$ from its ancestors to eliminate cross-terms.

2. Exact Marginalization (Theorem 3.2): The proof constructs the block matrix form of $(I - B)$ where $B$ is the structural matrix:

   \[I - B = \begin{pmatrix} I_S - B_{SS} & 0 \\ -B_{l,S} & 1 \end{pmatrix}\]

   Computing $(I-B)^T(I-B)$ and applying the Schur complement formula yields the exact marginalization property.

3. Expectation Gap Characterization (Proposition 3.4): The key inequality establishes:

   \[\Delta = \mathbb{E}[\text{Schur}(H(x))] - \text{Schur}(\mathbb{E}[H(x)]) = -\frac{1}{\sigma^2}\text{Cov}_{p(x)}(\nabla f_l(x))\]

   This uses the fact that sample-wise Hessians have constant diagonal $H_{ll}(x) = 1/\sigma^2$ for leaf nodes, while off-diagonal blocks contain the gradient $\nabla f_l(x)$.

4. Technical Innovation: The proof technique connecting Fisher information geometry to causal structure via Schur complements is novel. The localization property (Remark 3.2) showing that approximation errors concentrate only in parent neighborhoods is a sophisticated technical insight.

Datasets:

• Synthetic linear and non-linear ANMs with dimensions d ∈ [20, 1000]
• Sachs protein signaling network (d=12, N=5000)
• Sample sizes N ∈ [500, 5000]

Baselines: NOTEARS, DAGMA (continuous optimization), Recursive SCORE, Scalable SCORE (score-based methods)

Key Numbers:

• Linear ANMs: Zero edge violations across all dimensions (exact recovery)
• d=50 non-linear: SHD 2.6±2.4 vs DAGMA 93.0±12.6, 700× speedup (3.66s vs 2700s)
• d=100: SHD 64.6±17.3 vs Scalable SCORE 10.4±5.5, but 850× faster extraction (0.04s vs 34s)
• d=1000: 488.6±8.5 edge violations (approaching random baseline of 500), extraction time 0.56s
• Sachs dataset: SHD 13.8±1.2, TPR 0.40±0.04, 27× faster than DAGMA

Implementation: Neural score networks with Group Lasso regularization, Block-SSTS with tolerance γ=0.05, ridge regularization λ=10⁻⁴. Memory footprint O(d²), evaluated on single RTX 5090 GPU.

Limitations:

1. Additive Noise Assumption: RESTRICTIVE - Method fails completely under Post-Nonlinear (PNL) mechanisms where xi = (fi(xpa(i)) + εi)³, achieving only 3% TPR vs 95% under ANMs

2. Expectation Gap Accumulation: TECHNICAL - In non-linear ANMs, sequential marginalization accumulates approximation errors proportional to Cov(∇fl), degrading performance with graph depth

3. High-dimensional Statistical Limits: NATURAL - At d=1000 with N=5000, performance approaches random due to score network estimation variance in under-parameterized regimes

4. Homoscedastic Noise: TECHNICAL - Core theory assumes constant noise variance σ², though Appendix D.1 provides scale-invariant extensions

5. Dense Graph Performance: NATURAL - Method designed for sparse DAGs; performance degrades on dense graphs with many parallel structures

   Open Problems:

6. Extending to PNL Models: Developing exact Schur marginalization for arbitrary post-nonlinear transformations remains theoretically open

7. Finite-sample Theory: Establishing minimax optimal sample complexity bounds for SJIM estimation in high-dimensional regimes

1. Motivation: Spiking transformers achieve competitive accuracy with standard transformers while offering 38-57× energy efficiency on neuromorphic hardware. However, no theoretical framework exists to guide their design, creating a gap between empirical success and understanding of their computational limits.

2. Mathematical setup: Consider input sequences $X ∈ [0,1]^{n×d}$ with $n$ tokens and $d$ dimensions. Spiking neurons follow Leaky Integrate-and-Fire (LIF) dynamics:

   \[u_t = βu_{t-1} + I_t - v_{th}s_{t-1}\] \[s_t = Θ(u_t - v_{th})\]

   where $β ∈ (0,1)$ is membrane decay, $v_{th}$ is firing threshold, and $Θ(·)$ is Heaviside function. For spike rate encoding over $T$ timesteps with rate $x ∈ [0,1]$:

   \[E\left[\sum_{t=1}^T s_t\right] = xT\]

   Spiking self-attention processes spike trains $S^X ∈ {0,1}^{n×d×T}$ via:

   \[S^Q = SN(S^X W^Q), \quad S^K = SN(S^X W^K)\] \[A = \frac{1}{T}\sum_{t=1}^T S^Q_t (S^K_t)^⊤\]

   Assumptions:

   1. Compact domain $K ⊂ [0,1]^{n×d}$
   2. Lipschitz continuity with constant $L_f$
   3. Bounded spike rates $\bar{s} ∈ [ρ, 1-ρ]$ for some $ρ > 0$

3. Toy example: When $n=2, d=2$ and input $X = \begin{bmatrix} 0.8 & 0.2 \ 0.3 & 0.7 \end{bmatrix}$, spike rate encoding produces binary spike trains over $T$ timesteps where each dimension fires at rate proportional to its value.

4. Formal objective: Determine the minimum spike count $N_{total}$ required to achieve $ε$-approximation:

   \[\sup_{X ∈ K} \|f(X) - f_{spike}(X)\|_F < ε\]

The method provides explicit constructions for universal approximation via spike circuits:

1. Spike rate approximation (Lemma 1): For $x ∈ [ρ, 1-ρ]$ and error $δ$, require $T_0 = O(1/δ^2)$ timesteps

2. Exponential computation (Lemma 2): Compute $e^z$ via hierarchical spike coincidence using truncated Taylor series:

   \[e^z ≈ \sum_{j=0}^J \frac{z^j}{j!}\]

   with $j$-way coincidence circuits producing $s^{(C_j)}_t = \prod_{i=1}^j s^{(i)}_t$

3. Winner-Take-All circuit (Lemma 3): Implement softmax normalization $α_i = e_i/\sum_j e_j$ using lateral inhibition with global inhibitory feedback:

   \[I^{inh}(t) = \frac{1}{n}\sum_{j=1}^n s_j(t-1)\]

4. Universal construction (Theorem 7): Build spiking transformer with $L = O(\log(1/ε))$ layers, $H = O(n^2d/ε^2)$ heads, $T = O(n^2d_k^2/ε^2)$ timesteps

   Applying to toy example: For $2×2$ input, the WTA circuit creates 2 LIF neurons with cross-inhibition. Each neuron receives excitatory input proportional to $\exp(q^⊤k_i)$ and inhibitory feedback from the population sum, converging to normalized attention weights in $O(n^2/ε^2)$ timesteps.

Step 1 — Prior work:

• Yun et al. (2020) “Are transformers universal approximators of sequence-to-sequence functions?” proved standard transformers are universal approximators
• Maass (1997) “Networks of spiking neurons” established basic SNN expressivity for feedforward networks
• Zhou et al. (2023) “Spikformer” introduced practical spiking self-attention architectures

Step 2 — Delta: This paper extends universal approximation to the spiking domain with explicit circuit constructions, derives tight spike-count lower bounds via rate-distortion theory, and introduces input-dependent analysis using effective dimensions to explain the theory-practice gap.

Step 3 — Theory-specific assessment:

• Main theorem is predictable extension of existing universal approximation results to spiking domain
• Proof technique assembles known concentration bounds and rate-distortion theory in standard way
• Lower bound Ω(L_f^2 nd/ε^2) appears tight based on experiments, but no matching lower bound proof exists
• The effective dimension insight (d_eff ≪ nd) is the most novel contribution, explaining why practical systems need far fewer spikes

The circuit constructions are detailed but routine applications of spike encoding techniques. The rate-distortion derivation follows standard information theory.

Verdict: INCREMENTAL — solid theoretical foundation for existing empirical results, but techniques are straightforward extensions of known methods.

Main proof strategy uses four key components:

1. Concentration bounds for spike rate encoding: Apply Chernoff bound for bounded random variables:

   \[P[|\bar{s} - x| > δ] ≤ 2\exp(-2Tδ^2)\]

2. Hierarchical coincidence detection for exponentials: Encode $z$ as rate $r_z = (z+M)/(2M)$, compute powers via:

   \[E[\bar{s}^{(C_j)}] = r_z^j\]

   where $s^{(C_j)}_t = \prod_{i=1}^j s^{(i)}_t$

3. Winner-Take-All convergence analysis: Use Lyapunov function $V(r) = \sum_i (r_i - α_i)^2$ to prove exponential convergence of lateral inhibition dynamics

4. Rate-distortion lower bound derivation: For source $X \sim \text{Uniform}([0,1]^{nd})$ at distortion $D = (ε/L_f)^2$:

   \[R(D) = nd \cdot \log_2(L_f/ε)\]

   Combined with variance constraint requiring $T = Ω(L_f^2/ε^2)$ timesteps gives:

   \[E[N_{total}] = Ω(L_f^2 nd/ε^2)\]

5. Effective dimension reduction: Replace worst-case $nd$ with measured $d_{eff}$ based on PCA analysis capturing 95% variance, explaining practical efficiency.

   The key technical insight is using input distribution properties rather than worst-case bounds.

Datasets: CIFAR-10/100, ImageNet-1K (vision), SST-2 (NLP) Baselines: Spikformer, QKFormer, SpikingResformer, SpikingBERT Key numbers:

• Theoretical predictions validated with R² = 0.97 (p < 0.001)
• Observed/theoretical ratios: 2.08-2.63× across architectures
• Effective dimensions: deff = 47±3 (CIFAR-10), 68±4 (CIFAR-100), 89±7 (ImageNet)
• Design rule constant: C = 2.3 (95% CI: [1.9, 2.7])
• Energy efficiency: 38-57× vs standard transformers
• Timestep predictions within 2-18% error of optimal values

Implementation: SpikingJelly framework, AdamW optimizer, 10 random seeds for statistical rigor. All theoretical scaling laws validated experimentally.

Limitations:

1. Ideal spike rate encoding assumption - TECHNICAL (real systems use surrogate gradients with additional noise)
2. WTA circuit requires non-zero input separation - TECHNICAL (can be relaxed with regularization)
3. Compact domain restriction [0,1]^{n×d} - NATURAL (standard in universal approximation theory)
4. LIF neuron model only - TECHNICAL (other spiking models could be analyzed similarly)

5. Static effective dimension measurement - RESTRICTIVE (deff may vary across samples/tasks)

   Open problems:

6. Extend analysis to temporal coding schemes beyond rate encoding (e.g., STDP-based learning)
7. Characterize spike-count requirements for sparse attention patterns emerging from lateral inhibition dynamics