Spectral Bias in Neural Networks and Resolution via Fourier Kernel Methods

Table of Contents

• 1. Introduction
• 2. The Spectral Bias Phenomenon
  • 2.1 Empirical Evidence
  • 2.2 Formal Definition and the Frequency Principle
  • 2.3 Fourier Spectrum of ReLU Networks
• 3. NTK Analysis of Spectral Bias
  • 3.1 Neural Tangent Kernel as Kernel Regression
  • 3.2 Eigendecomposition and Frequency-Dependent Convergence
  • 3.3 Spherical Harmonics and the Decay of NTK Eigenvalues
• 4. Fourier Feature Maps
  • 4.1 Bochner’s Theorem
  • 4.2 Rahimi and Recht: Random Fourier Features
  • 4.3 Positional Encodings as Deterministic Fourier Features
• 5. Effect on the NTK Spectrum
  • 5.1 The Composed Kernel
  • 5.2 Stationarity and Bandwidth Control
  • 5.3 Eigenvalue Equalization and Kernel Preconditioning
• 6. Convergence Rate Analysis
  • 6.1 Baseline Convergence in the NTK Regime
  • 6.2 Convergence After Fourier Feature Lifting
• 7. Applications
  • 7.1 NeRF Positional Encodings
  • 7.2 SIREN: Sinusoidal Activation Networks
  • 7.3 Occupancy Networks and Other Implicit Representations
• 8. Cross-Attention as Adaptive Frequency Selection
  • 8.1 The Limitation of Fixed Fourier Features
  • 8.2 Multiscale Token Bank
  • 8.3 Cross-Attention over Frequency Tokens
  • 8.4 Adaptive Frequency Masking
  • 8.5 Adaptive Frequency Enhancement via DFT Injection
  • 8.6 Dual-Network Architecture for PDEs
  • 8.7 Experimental Results
• 8.8 PyTorch Implementation
• 9. References

1. Introduction 🗺️

Implicit neural representations (INRs) are networks that map coordinates \((x, y, z, t, \ldots)\) directly to signal values: pixel intensities, signed distances, density fields, or radiance. The appeal is compact parameterization — a single MLP encodes a whole scene or object.

The obstacle is fundamental: plain MLPs with smooth (e.g., ReLU, tanh) activations exhibit a strong spectral bias — they learn to represent low-frequency content orders of magnitude faster than high-frequency content. For natural signals (images, 3D geometry, audio), the high-frequency components carry texture, edges, and fine detail. A spectrally biased network trained on a 3D scene will render a blurry blob long before it can reproduce surface texture.

This note develops the spectral bias story from three angles:

1. Empirical and Fourier-analytic (Rahaman et al., 2019): what we observe in practice and what the Fourier transform of ReLU networks reveals.
2. NTK-theoretic (Jacot et al., 2018; Basri et al., 2019): how the neural tangent kernel’s eigenspectrum directly encodes the frequency-learning schedule.
3. Feature-map resolution (Rahimi & Recht, 2007; Tancik et al., 2020): how Fourier feature maps reshape the induced kernel, equalizing eigenvalues across frequency modes and recovering high-frequency learning ability.

2. The Spectral Bias Phenomenon 📊

2.1 Empirical Evidence

The following experiment (following Rahaman et al., 2019) cleanly isolates the phenomenon.

Setup. Train a deep ReLU network \(f_\theta : \mathbb{R} \to \mathbb{R}\) to regress a target \(g(x) = \sum_{k} A_k \sin(2\pi k x)\) on \(x \in [0,1]\), where frequencies \(k \in \{1, 5, 10, 40\}\) are simultaneously present.

Observation. Monitoring the residual power spectrum \(|\hat{r}(k)|^2\) where \(\hat{r}(k) = \widehat{g - f_{\theta^{(t)}}}(k)\) at training iteration \(t\):

• The \(k=1\) component vanishes within a few hundred iterations.
• The \(k=5\) component decays an order of magnitude later.
• \(k=10, 40\) components remain present long into training.

Crucially, this ordering holds even when higher-frequency amplitudes \(A_k\) are larger. The network’s preference is not amplitude-driven but frequency-driven.

A secondary experiment perturbs the parameters of a trained network by additive Gaussian noise \(\epsilon \sim \mathcal{N}(0, \sigma^2 I)\) and measures the effect on each frequency component. Low-frequency components are nearly unaffected; high-frequency components degrade rapidly. This shows that high-frequency representations occupy a small volume in parameter space — they require precise parameter coordination, making them both hard to learn and fragile.

2.2 Formal Definition and the Frequency Principle

Definition (Spectral Bias). Let \(f^{(t)} : \mathbb{R}^d \to \mathbb{R}\) be the network function at training time \(t\), and let \(f^* : \mathbb{R}^d \to \mathbb{R}\) be the target function. The network exhibits spectral bias if the residual error \(r^{(t)} = f^* - f^{(t)}\) has the property that, for any \(\xi_1, \xi_2 \in \mathbb{R}^d\) with \(\|\xi_1\| < \|\xi_2\|\):

\[|\hat{r}^{(t)}(\xi_2)| \gg |\hat{r}^{(t)}(\xi_1)|\]

throughout most of training, i.e., the high-frequency residual decays far slower than the low-frequency residual.

Definition (Frequency Principle / F-Principle). Training dynamics satisfy the frequency principle if the Fourier components of the training residual decay approximately in order of increasing frequency \(\|\xi\|\).

The two definitions are closely related. The frequency principle is the dynamic statement (frequencies learned in order); spectral bias is the static description of the outcome (networks are biased toward low-frequency solutions).

2.3 Fourier Spectrum of ReLU Networks

Rahaman et al. (2019) provide a Fourier-analytic characterization of what ReLU networks can express and at what cost.

Setup. Let \(f_\theta : \mathbb{R}^d \to \mathbb{R}\) be an \(L\)-layer ReLU network with width \(n\) and weight matrices \(\{W^{(l)}\}\). A ReLU network is piecewise linear, partitioning its input domain into a set of linear regions (convex polytopes) \(\{P_i\}\), on each of which \(f_\theta\) is an affine function.

Theorem 1 (Rahaman et al., 2019 — Fourier Spectrum of ReLU Networks). For a piecewise-linear function \(f_\theta\) with linear regions \(\{P_i\}\), the Fourier transform satisfies:

\[\hat{f}_\theta(\xi) = \sum_i \sum_{n=0}^{d} \frac{C_n(\theta, \xi) \cdot \mathbf{1}_{H^i_n}(\xi)}{\|\xi\|^{n+1}}\]

where \(H^i_n\) is the set of frequency directions aligned with \(n\)-dimensional faces of polytope \(P_i\), and \(C_n(\theta, \xi)\) is a bounded coefficient depending on the local affine map and face geometry.

The upshot: the Fourier transform of a ReLU network decays as \(\|\xi\|^{-(n+1)}\) for some \(n \in \{0, \ldots, d\}\), so the spectral content is suppressed polynomially at high frequencies. For a generic direction, \(n = 0\) and the decay is \(\|\xi\|^{-1}\) — the mildest suppression possible. But for directions perpendicular to high-dimensional faces, \(n = d-1\) and the decay is \(\|\xi\|^{-d}\), much stronger.

Corollary. The rate at which the residual’s \(\xi\)-component changes satisfies:

\[\left|\frac{d\hat{r}^{(t)}(\xi)}{dt}\right| = O\!\left(\|\xi\|^{-\Delta}\right), \quad \Delta \in [1, d]\]

meaning high-frequency residual components change more slowly during gradient descent, giving the formal basis for the empirical observations in §2.1.

3. NTK Analysis of Spectral Bias 📐

3.1 Neural Tangent Kernel as Kernel Regression

The neural tangent kernel (NTK), introduced by Jacot et al. (2018), characterizes the infinite-width limit of neural network training dynamics.

Definition (NTK). For a network \(f(\cdot;\theta) : \mathbb{R}^d \to \mathbb{R}\) with parameters \(\theta \in \mathbb{R}^P\), the neural tangent kernel is:

\[K_{\mathrm{NTK}}(x, x') = \left\langle \frac{\partial f(x;\theta)}{\partial \theta},\, \frac{\partial f(x';\theta)}{\partial \theta} \right\rangle_{\mathbb{R}^P}\]

In the infinite-width limit under gradient flow on a squared loss \(\mathcal{L} = \frac{1}{2}\|f_\theta(X) - y\|^2\), two facts hold:

1. Constancy: \(K_{\mathrm{NTK}}\) remains constant throughout training (i.e., \(\frac{d}{dt} K_{\mathrm{NTK}} = 0\)).
2. Kernel regression limit: The network predictions \(\hat{y}^{(t)} = f_\theta^{(t)}(X_{\mathrm{test}})\) converge to the kernel regression solution:

\[\hat{y}^{(t)}_{\mathrm{test}} = K_{\mathrm{test}} K^{-1} \left(I - e^{-\eta K t}\right) y\]

where \(K = K_{\mathrm{NTK}}(X_{\mathrm{train}}, X_{\mathrm{train}}) \in \mathbb{R}^{n \times n}\) is the train–train kernel matrix, \(K_{\mathrm{test}} = K_{\mathrm{NTK}}(X_{\mathrm{test}}, X_{\mathrm{train}}) \in \mathbb{R}^{m \times n}\), and \(\eta\) is the learning rate.

3.2 Eigendecomposition and Frequency-Dependent Convergence

Let \(K = Q \Lambda Q^\top\) be the eigendecomposition with \(\Lambda = \mathrm{diag}(\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_n > 0)\) and orthonormal eigenvectors \(Q = [q_1 | q_2 | \cdots | q_n]\).

Proposition (NTK Convergence in Eigenbasis). Under gradient flow, the error in the eigenbasis satisfies:

\[Q^\top \!\left(\hat{y}^{(t)}_{\mathrm{train}} - y\right) = -e^{-\eta \Lambda t}\, Q^\top y\]

That is, the \(i\)-th eigenvector component of the error decays as \(e^{-\eta \lambda_i t}\).

Proof sketch. The training residual \(r^{(t)} = y - \hat{y}^{(t)}_{\mathrm{train}}\) evolves as \(\dot{r}^{(t)} = -K r^{(t)}\), giving \(r^{(t)} = e^{-Kt} r^{(0)}\). Rotating by \(Q^\top\): \((Q^\top r^{(t)})_i = e^{-\lambda_i t} (Q^\top r^{(0)})_i\). \(\square\)

Definition (Convergence Rate of Mode \(i\)). The \(i\)-th mode converges at rate \(\eta \lambda_i\); its half-life is \((\eta \lambda_i)^{-1} \log 2\).

This is the core of the NTK-based explanation of spectral bias:

💡 If \(\lambda_i \ll \lambda_j\) for eigenvectors \(q_i, q_j\) corresponding to high and low frequencies respectively, then the high-frequency mode takes far longer to fit.

Key question: Do the NTK eigenvectors actually align with frequency modes, and do low-frequency modes have larger eigenvalues?

3.3 Spherical Harmonics and the Decay of NTK Eigenvalues

For training data distributed uniformly on the \((d-1)\)-sphere \(\mathbb{S}^{d-1}\), the NTK matrix forms a zonal kernel — it depends only on the inner product \(x^\top x'\):

\[K_{\mathrm{NTK}}(x, x') = h(x^\top x'), \quad \|x\| = \|x'\| = 1\]

By the theory of zonal kernels on spheres, the eigenfunctions of the associated integral operator

\[(\mathcal{K} f)(x) = \int_{\mathbb{S}^{d-1}} h(x^\top x') f(x')\, d\sigma(x')\]

are exactly the spherical harmonics \(\{Y_{l,j}\}\) of degree \(l = 0, 1, 2, \ldots\) and order \(j = 1, \ldots, N(d,l)\), where \(N(d,l)\) is the multiplicity of degree \(l\).

Theorem (Basri et al., 2019; Bietti & Mairal, 2019). For a two-layer overparameterized ReLU network with weights initialized from a standard isotropic Gaussian, the NTK eigenvalue \(\mu_l\) corresponding to degree-\(l\) spherical harmonics satisfies:

\[\mu_l \sim l^{-(d+1)}\]

as \(l \to \infty\). That is, eigenvalues decay polynomially in frequency, with exponent \(d+1\).

Corollary (Frequency-Dependent Convergence Rate). In the NTK regime with uniform data on \(\mathbb{S}^{d-1}\), the degree-\(l\) spherical harmonic component of the training error decays as \(e^{-\eta \mu_l t} \approx e^{-\eta l^{-(d+1)} t}\). For the error to reach half its initial value requires time:

\[t_{1/2}(l) \sim l^{d+1}\]

This is the formal quantification of spectral bias: high-frequency modes (large \(l\)) take polynomially longer to learn, with the gap growing with input dimension \(d\).

4. Fourier Feature Maps 🔑

4.1 Bochner’s Theorem

The theoretical bridge between shift-invariant kernels and random feature maps is Bochner’s theorem from harmonic analysis.

Definition (Shift-Invariant Kernel). A kernel \(k : \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}\) is shift-invariant (stationary) if \(k(x, x') = k(x - x')\) for some function \(k : \mathbb{R}^d \to \mathbb{R}\).

Theorem (Bochner). A continuous, bounded, shift-invariant kernel \(k : \mathbb{R}^d \to \mathbb{R}\) is positive definite if and only if it is the Fourier transform of a finite non-negative Borel measure \(p\) on \(\mathbb{R}^d\):

\[k(x - x') = \int_{\mathbb{R}^d} p(\omega)\, e^{i\omega^\top(x - x')}\, d\omega\]

If \(k(0) = 1\) (normalized), then \(p\) is a probability density on \(\mathbb{R}^d\) and is called the spectral density of the kernel.

Proof sketch. The forward direction (positive definite \(\Rightarrow\) Fourier transform of measure) follows from the Herglotz–Bochner theorem in harmonic analysis. The converse is the Bochner–Schwartz theorem: if \(p \geq 0\), then \(k(x - x') = \mathbb{E}_\omega[e^{i\omega^\top(x-x')}]\) is a convex combination of positive definite functions \(\omega \mapsto e^{i\omega^\top \cdot}\), hence positive definite. \(\square\)

Corollary (Real-Valued Form). Since \(k\) is real and symmetric, the imaginary parts cancel in the integral:

\[k(x - x') = \int_{\mathbb{R}^d} p(\omega) \cos(\omega^\top(x - x'))\, d\omega = \mathbb{E}_{\omega \sim p}\!\left[\cos(\omega^\top(x - x'))\right]\]

This is the key identity: the kernel is an expectation over cosine features.

4.2 Rahimi and Recht: Random Fourier Features

Rahimi & Recht (2007) exploit Bochner’s theorem to construct a randomized, finite-dimensional feature map that approximates any shift-invariant kernel.

Definition (Random Fourier Feature Map). Given a shift-invariant kernel \(k\) with spectral density \(p\), sample \(\omega_1, \ldots, \omega_m \overset{\mathrm{iid}}{\sim} p\) and define the random Fourier feature map \(z : \mathbb{R}^d \to \mathbb{R}^{2m}\):

\[z(x) = \frac{1}{\sqrt{m}}\!\left[\cos(\omega_1^\top x),\, \sin(\omega_1^\top x),\, \ldots,\, \cos(\omega_m^\top x),\, \sin(\omega_m^\top x)\right]^\top\]

Proposition (Unbiased Kernel Approximation). The inner product \(z(x)^\top z(x')\) is an unbiased estimator of the kernel:

\[\mathbb{E}\!\left[z(x)^\top z(x')\right] = k(x - x')\]

Proof. Using the Euler identity \(\cos(\alpha)\cos(\beta) + \sin(\alpha)\sin(\beta) = \cos(\alpha - \beta)\):

\[\mathbb{E}\!\left[z(x)^\top z(x')\right] = \mathbb{E}_{\omega \sim p}\!\left[\cos(\omega^\top x)\cos(\omega^\top x') + \sin(\omega^\top x)\sin(\omega^\top x')\right] = \mathbb{E}_\omega\!\left[\cos(\omega^\top(x - x'))\right] = k(x - x') \quad \square\]

Concentration. By Hoeffding’s inequality applied to the bounded random variables (since \(|\cos(\cdot)| \leq 1\)), the approximation error concentrates:

\[\Pr\!\left[\sup_{x, x'} \left|z(x)^\top z(x') - k(x - x')\right| > \epsilon\right] \leq C \cdot e^{-m\epsilon^2/C'}\]

so \(m = O(\epsilon^{-2} \log(1/\delta))\) features suffice for uniform \(\epsilon\)-approximation with probability \(1-\delta\).

Example: RBF Kernel. For the Gaussian (RBF) kernel \(k(x - x') = \exp(-\|x-x'\|^2 / 2\sigma^2)\), the spectral density is \(p(\omega) = \mathcal{N}(0, \sigma^{-2} I_d)\). Sampling \(\omega \sim \mathcal{N}(0, \sigma^{-2} I)\) and forming \(z\) gives an unbiased approximation to the RBF kernel.

4.3 Positional Encodings as Deterministic Fourier Features

Tancik et al. (2020) make the connection between Fourier features and neural network conditioning explicit and replace Monte Carlo sampling with a fixed, deterministic feature map.

Definition (Fourier Feature Map / Positional Encoding). Given a frequency matrix \(B \in \mathbb{R}^{m \times d}\) with rows \(b_1, \ldots, b_m\), define:

\[\gamma(v) = \left[\cos(2\pi b_1^\top v),\, \sin(2\pi b_1^\top v),\, \ldots,\, \cos(2\pi b_m^\top v),\, \sin(2\pi b_m^\top v)\right]^\top \in \mathbb{R}^{2m}\]

The MLP then operates on \(\gamma(v)\) rather than \(v\) directly: \(f_\theta(\gamma(v))\).

Three variants arise depending on how \(B\) is chosen:

The NeRF encoding maps \(v \mapsto [\cos(v), \sin(v), \cos(2v), \sin(2v), \ldots, \cos(2^{L-1}v), \sin(2^{L-1}v)]\) — a deterministic frequency ladder spanning from fundamental to \(2^{L-1}\) times the fundamental.

5. Effect on the NTK Spectrum 🔬

5.1 The Composed Kernel

When inputs are lifted via \(\gamma\), the effective NTK of the network \(f_\theta(\gamma(\cdot))\) changes. Let \(h_{\mathrm{NTK}} : \mathbb{R} \to \mathbb{R}\) denote the NTK of the base MLP (which, for typical architectures, is a function of the inner product \(\gamma(v)^\top \gamma(v')\)).

Definition (Composed Kernel). The composed kernel induced by the Fourier feature map is:

\[K_\gamma(v, v') = h_{\mathrm{NTK}}\!\left(\gamma(v)^\top \gamma(v')\right) = h_{\mathrm{NTK}}\!\left(\sum_{j=1}^m \cos(b_j^\top(v - v'))\right)\]

Since the argument \(\gamma(v)^\top \gamma(v') = h_\gamma(v - v')\) depends only on \(v - v'\), the composed kernel \(K_\gamma\) is shift-invariant (stationary).

Theorem (Tancik et al., 2020 — Stationarity). If \(h_{\mathrm{NTK}}\) is a function of the inner product and \(\gamma\) is a Fourier feature map, then the composed kernel \(K_\gamma(v, v') = h_{\mathrm{NTK}}(h_\gamma(v - v'))\) is shift-invariant. As a composition of an isotropic function with a shift-invariant kernel, it is itself shift-invariant, and by Bochner’s theorem has a spectral density \(\tilde{p}\) satisfying:

\[K_\gamma(v - v') = \int_{\mathbb{R}^d} \tilde{p}(\omega) e^{i\omega^\top(v - v')}\, d\omega\]

5.2 Stationarity and Bandwidth Control

The original NTK of a coordinate MLP is not stationary — it decays as inputs move apart, with a decay rate that depends strongly on frequency. The feature map \(\gamma\) converts it into a stationary kernel whose bandwidth is governed by the scale of \(B\).

Key mechanism: For the Gaussian random variant with \(b_j \sim \mathcal{N}(0, \sigma^2 I)\), the feature map inner product approximates a Gaussian kernel of bandwidth \(\propto \sigma^{-1}\):

\[h_\gamma(v - v') \approx \exp\!\left(-\frac{\sigma^2 \|v - v'\|^2}{2}\right)\]

Composing with \(h_{\mathrm{NTK}}\) (which is an increasing function near \(1\) for typical normalized inputs), the composed kernel has a width set by \(\sigma\). Larger \(\sigma\) makes the feature map project to higher frequencies and narrows the kernel in spatial domain, corresponding to a kernel that resolves finer detail.

The training dynamics (§3.1) with the composed kernel are:

\[\hat{y}^{(t)}_{\mathrm{test}} = K^\gamma_{\mathrm{test}} (K^\gamma)^{-1}(I - e^{-\eta K^\gamma t}) y\]

where \(K^\gamma\) is the \(n \times n\) composed kernel matrix on training points.

5.3 Eigenvalue Equalization and Kernel Preconditioning

The central claim of Tancik et al. (2020) is that the Fourier feature map equalizes the eigenvalues of the induced kernel across frequency modes, allowing all modes to converge at comparable rates.

Heuristic argument for equalization. The original coordinate MLP NTK has eigenvalues decaying as \(\mu_l \sim l^{-(d+1)}\) for spherical harmonic degree \(l\) (§3.3). After lifting via \(\gamma\), the composed kernel \(K_\gamma\) is a stationary kernel whose spectral density \(\tilde{p}\) is shaped by the distribution of frequencies \(b_j\).

For the random Fourier feature map with \(b_j \sim \mathcal{N}(0, \sigma^2 I)\):

• The spectral density \(\tilde{p}\) is supported over a ball of radius \(\sim \sigma\) in frequency space.
• Modes with \(\|\xi\| \lesssim \sigma\) all receive comparable spectral weight.
• The eigenvalue decay, instead of being \(l^{-(d+1)}\), is approximately flat for \(l \lesssim \sigma\) and drops sharply for \(l \gg \sigma\).

This is the sense in which the Fourier feature map preconditions the kernel: it converts the strongly non-uniform eigenvalue distribution of the raw NTK into a more uniform one over the frequency range \([0, \sigma]\).

Connection to kernel preconditioning. Preconditioned gradient descent on a kernel regression problem with kernel \(K\) and preconditioner \(P\) is equivalent to vanilla gradient descent with kernel \(P^{1/2} K P^{1/2}\). The Fourier feature map effectively implements a preconditioning that boosts high-frequency eigenvalues, making the optimization landscape more isotropic.

6. Convergence Rate Analysis 📈

6.1 Baseline Convergence in the NTK Regime

Recall the NTK training dynamics (§3.1): the \(i\)-th eigenvector component of the training residual decays as

\[[Q^\top r^{(t)}]_i = e^{-\eta \lambda_i t} [Q^\top r^{(0)}]_i\]

Definition (Frequency Mode Convergence Rate). For a target \(f^*\) with dominant Fourier mode \(\omega\), let \(\lambda_\omega\) be the NTK eigenvalue associated with the eigenvector most aligned with frequency \(\omega\). The convergence rate of mode \(\omega\) is \(\eta \lambda_\omega\), and the mode fits to accuracy \(\epsilon\) at time:

\[t_\epsilon(\omega) = \frac{\log(1/\epsilon)}{\eta \lambda_\omega}\]

For a plain coordinate MLP, \(\lambda_\omega \sim \|\omega\|^{-(d+1)}\), so:

\[\boxed{t_\epsilon(\omega) \sim \frac{\|\omega\|^{d+1} \log(1/\epsilon)}{\eta}}\]

The convergence time grows as the \((d+1)\)-th power of frequency. In 3D (\(d=3\)), doubling the spatial frequency increases the required training time by \(2^4 = 16\times\).

6.2 Convergence After Fourier Feature Lifting

After lifting \(v \mapsto \gamma(v)\) with \(B \sim \mathcal{N}(0, \sigma^2 I)\), the composed kernel \(K_\gamma\) has an approximately flat spectrum for \(\|\omega\| \lesssim \sigma\). In this regime, all eigenvalues satisfy \(\tilde\lambda_\omega \approx c\) for some constant \(c\) (the plateau value of the spectral density), and:

\[t_\epsilon(\omega) \approx \frac{\log(1/\epsilon)}{\eta c}, \quad \|\omega\| \lesssim \sigma\]

Convergence is frequency-independent for all modes within the bandwidth.

Theorem (Tancik et al., 2020 — informal). For an MLP trained on coordinates mapped through \(\gamma(v) = [\cos(Bv), \sin(Bv)]\) with \(B \sim \mathcal{N}(0, \sigma^2 I)\):

1. The composed NTK is stationary with bandwidth parameter \(\sigma\).
2. For target frequencies \(\|\omega\| \leq \sigma\), the convergence rate is \(\Theta(\tilde\lambda_\text{max})\), independent of \(\|\omega\|\).
3. For target frequencies \(\|\omega\| > \sigma\), the mode lies outside the kernel’s effective support and cannot be learned.

Quantitative comparison. Consider frequency \(\|\omega\| = L\) in \(d\) dimensions:

The Fourier feature map converts polynomial-in-frequency convergence time to constant convergence time, for all modes within the kernel bandwidth.

7. Applications 🖼️

7.1 NeRF Positional Encodings

Neural Radiance Fields (NeRF; Mildenhall et al., 2020) represent a 3D scene as a continuous volumetric function \(F_\theta : (\mathbf{x}, \mathbf{d}) \to (\mathbf{c}, \sigma)\) mapping 3D coordinates \(\mathbf{x} \in \mathbb{R}^3\) and viewing direction \(\mathbf{d} \in \mathbb{S}^2\) to RGB color \(\mathbf{c}\) and volume density \(\sigma\).

Without positional encoding, an MLP trained on \((\mathbf{x}, \mathbf{d})\) directly suffers spectral bias and produces blurry, low-frequency renderings. The solution (motivated by, and empirically validated prior to, the Tancik et al. theory) is the NeRF positional encoding:

\[\gamma(p) = \left(\sin(2^0 \pi p),\, \cos(2^0 \pi p),\, \sin(2^1 \pi p),\, \cos(2^1 \pi p),\, \ldots,\, \sin(2^{L-1} \pi p),\, \cos(2^{L-1} \pi p)\right)\]

applied component-wise to each coordinate. With \(L = 10\) for position and \(L = 4\) for direction, NeRF can recover fine textures, specular highlights, and thin geometric details that a plain MLP cannot.

The Tancik et al. analysis provides the retrospective explanation: the log-spaced encoding ensures the composed kernel has substantial spectral support out to frequency \(2^{L-1}\), with the geometric progression providing roughly uniform (in log-frequency) coverage.

7.2 SIREN: Sinusoidal Activation Networks

A distinct approach to the spectral bias problem is SIREN (Sitzmann et al., 2020): replace all ReLU activations with \(\sin(\omega_0 \cdot)\) for a hyperparameter \(\omega_0 > 0\).

Mechanism. A sine-activation network computes:

\[f_\theta(x) = W_L \sin(\omega_0 (W_{L-1} \sin(\omega_0 (\cdots W_1 \sin(\omega_0 W_0 x + b_0) \cdots)) + b_{L-1})) + b_L\]

Unlike piecewise-linear networks, sine activations are infinitely differentiable and have Fourier transforms supported everywhere — \(\hat{\sigma}(\xi) = \delta(\xi - \omega_0) / 2 + \delta(\xi + \omega_0) / 2\). Compositions create superpositions of sinusoids at frequencies \(\omega_0^k\) for \(k = 1, \ldots, L\), enabling the network to represent high-frequency functions without needing fine-tuned weight cooperation.

Surprisingly, the Fourier-theoretic argument is complementary to the Fourier features argument: SIREN changes the architecture to natively produce high-frequency outputs, whereas Fourier features change the input representation to reshape the NTK.

Practical considerations. SIREN requires a careful initialization scheme (weights initialized as \(\mathcal{U}(-\sqrt{6/n}, \sqrt{6/n})\) with \(\omega_0\) chosen to control frequency support) and is sensitive to these choices. Fourier feature maps, by contrast, can be bolted onto any existing MLP architecture.

7.3 Occupancy Networks and Other Implicit Representations

Occupancy networks (Mescheder et al., 2019) represent 3D geometry as \(f_\theta : \mathbb{R}^3 \to [0,1]\) predicting whether a point is inside the object. Without positional encoding, occupancy networks learn coarse geometry (large-scale shape) but miss fine surface details (cracks, ridges, engravings).

Fourier feature maps (or learned embeddings) applied to the input coordinates of occupancy networks enable representation of fine-grained geometry — thin structures, sharp creases — by elevating the NTK’s effective bandwidth.

Empirical comparison across implicit representation tasks (Tancik et al., 2020):

8. Cross-Attention as Adaptive Frequency Selection 🎯

8.1 The Limitation of Fixed Fourier Features

The Tancik et al. approach (§4–§6) selects the frequency matrix \(B\) once at initialization — the bandwidth \(\sigma\) is a global hyperparameter fixed before training begins. This raises a fundamental question: what if the target function has spatially non-uniform frequency content? A texture with fine detail in one region and smooth gradients in another would ideally use high \(\sigma\) locally and low \(\sigma\) elsewhere, but a fixed \(B\) cannot adapt.

Moreover, even for globally high-frequency targets, random sampling of \(B \sim \mathcal{N}(0, \sigma^2 I)\) may wastefully duplicate nearby frequency slots while missing informative frequency bands. The random feature bank is input-independent — every query to the network uses the same frequency dictionary regardless of where in the domain we are evaluating.

The proposal of Feng et al. (2025): treat the Fourier feature bank as a key–value memory and use cross-attention to let the network query it adaptively at every layer. The attention weights become learned, input-dependent frequency selectors, dynamically routing gradient signal to the frequency components the current residual requires.

8.2 Multiscale Token Bank 🏗️

Definition (Multiscale Random Fourier Feature Bank). Sample base frequencies \(\omega_m \sim \mathcal{N}(0, \sigma^{-2} I_{d_\text{in}})\) for \(m = 1, \ldots, M_\text{base}\). Form dyadic dilations:

\[\tilde{\omega}_{m,k} = 2^k \omega_m, \quad k = 0, 1, \ldots, K\]

giving \(M = M_\text{base}(K+1)\) total frequencies spanning \(K+1\) octaves. Attach learnable amplitude envelopes:

\[a_{m,k} = \exp(-\beta \|\tilde{\omega}_{m,k}\|_2), \quad \beta \geq 0 \text{ (trainable)}\]

so the network can suppress or amplify each octave. The full feature vector is:

\[\phi(x) = \sqrt{\frac{1}{M}} \bigl[a_{m,k} \cos(\tilde{\omega}_{m,k}^\top x + b_{m,k})\bigr]_{(m,k)} \in \mathbb{R}^M\]

Token bank. Reshape \(\phi(x)\) with token width \(d_q\) to obtain the frequency token bank:

\[H(x) \in \mathbb{R}^{N_\text{tok} \times d_q}, \quad N_\text{tok} = M / d_q\]

Each row \(H_j(x)\) is a \(d_q\)-dimensional token encoding a group of adjacent frequency channels evaluated at \(x\). Crucially, \(H(x)\) is input-dependent — the cosine evaluations change with \(x\), so different inputs present different token activations to the attention mechanism.

8.3 Cross-Attention over Frequency Tokens 🔑

The main network is a stack of \(L\) query layers \(Q^{(0)}(x), \ldots, Q^{(L)}(x) \in \mathbb{R}^{d_q}\). At each layer \(l\), the query attends over the token bank \(H(x)\):

\[Q_l = Q^{(l)}(x) W_Q^{(l)}, \quad K_l = H(x) W_K^{(l)}, \quad V_l = H(x) W_V^{(l)}\]

with \(W_Q^{(l)}, W_K^{(l)}, W_V^{(l)} \in \mathbb{R}^{d_q \times d_q}\). The cross-attention output is:

\[\mathrm{CA}(Q^{(l)}(x),\, H(x)) = \mathrm{softmax}\!\left(\frac{Q_l K_l^\top}{\sqrt{d_q}}\right) V_l \in \mathbb{R}^{d_q}\]

This is then added residually, followed by a feedforward update:

\[\tilde{Q}^{(l)}(x) = Q^{(l)}(x) + \mathrm{CA}(Q^{(l)}(x),\, H(x))\]

\[Q^{(l+1)}(x) = \tilde{Q}^{(l)}(x) + \sigma\!\left(W^{(l)} \tilde{Q}^{(l)}(x) + b^{(l)}\right)\]

Why does this help? In the language of §3–§5, the fixed Fourier feature map \(\gamma(x)\) flattens the NTK eigenspectrum uniformly (§5.3). The cross-attention weights \(\mathrm{softmax}(Q_l K_l^\top / \sqrt{d_q})\) form a learned, input-dependent reweighting of the frequency tokens. At layer \(l\), the network concentrates attention mass on the frequency tokens \(H_j(x)\) that are most useful for reducing the current residual — effectively performing online eigenvalue reweighting beyond what the static feature map achieves.

8.4 Adaptive Frequency Masking 🎭

Early in training, the model has not yet learned which frequency tokens are informative; random attention over all \(N_\text{tok}\) tokens adds noise. Feng et al. introduce adaptive frequency masking: augment the attention logits with a learnable mask \(M^{(l)}\):

\[A^{(l)} = \frac{Q_l K_l^\top}{\sqrt{d_q}} + M^{(l)}, \quad M^{(l)} = [0;\; \eta_l \cdot \mathbf{1}]\]

where \(\eta_l \leq 0\) is initialized to a large negative value (suppressing the “posterior” tokens — those encoding the finer-scale octaves) and is gradually relaxed toward \(0\) during training.

Effect: the model begins by attending only to the coarse-frequency tokens (spectral bias is still present but at the token level), then progressively unlocks finer-scale tokens as training proceeds. This implements a curriculum over frequency scales, preventing the optimizer from wasting early gradient steps on high-frequency tokens that the coarse network cannot yet exploit.

8.5 Adaptive Frequency Enhancement via DFT Injection 📡

Fixed random features can miss sharp, localized frequency peaks in the target. The Adaptive Frequency Enhancement (AFE) module detects such peaks during training and injects them directly.

Procedure. At a designated checkpoint, compute the Discrete Fourier Transform of the current network output (or residual) on a reference grid:

\[\hat{u}_{\theta,k} = \mathrm{DFT}[u_\theta]_k\]

Identify the dominant unresolved modes:

\[\zeta = \max_{k \in \mathcal{B}} |\hat{u}_{\theta,k}|, \quad \mathcal{K}_\text{post} = \{k \in \mathcal{B} : |\hat{u}_{\theta,k}| > \lambda \zeta\}, \quad 0 < \lambda < 1\]

Construct posterior tokens exactly tuned to these modes:

\[\phi_\text{post}(x) = \sqrt{\frac{2}{M_\text{post}}} \cos(\Omega_\text{post} x + b_\text{post}), \quad \omega_k^\text{post} = 2k\pi\]

and augment the token bank:

\[H_\text{aug}(x) = [H_\text{base}(x);\; H_\text{post}(x)] \in \mathbb{R}^{(n_\text{base} + n_\text{post}) \times d_q}\]

No architectural modification is needed — the augmented bank slots into the cross-attention mechanism directly, since the attention is over the row dimension of \(H\). The newly injected tokens carry exact frequency information that random sampling may have missed.

8.6 Dual-Network Architecture for PDEs 📐

For physics-informed neural networks (PINNs), the spectral bias manifests as slow convergence of the high-frequency PDE solution components. Feng et al. propose decomposing the solution as:

\[u(x;\theta) = u_h(x;\theta_h) + \alpha\, u_l(x;\theta_l)\]

where \(u_h\) (with the full cross-attention architecture) captures high-frequency components and \(u_l\) (a lightweight standard MLP) captures the smooth baseline. The scalar \(\alpha\) is a trainable mixing factor.

Training loss. For a PDE \(\mathcal{N}[u] = f\) on \(\Omega\) with boundary condition \(u = g\) on \(\partial\Omega\):

\[\mathcal{L} = \int_\Omega \bigl(\mathcal{N}[u_h + \alpha u_l](x) - f(x)\bigr)^2 \rho_r(x)\,dx + \gamma \int_{\partial\Omega} \bigl((u_h - g)^2 + \alpha^2 u_l^2\bigr) \rho_b(x)\,dx\]

with \(u_l\) satisfying a homogeneous BC \(u_l|_{\partial\Omega} = 0\) so boundary conditions are handled by \(u_h\) alone.

Optimal scaling. For linear operators \(\mathcal{N}\), the optimal \(\alpha\) admits a closed form:

\[\alpha_\text{opt} = -\frac{\int_\Omega \bigl(\mathcal{N}[u_h] - f\bigr)\,\mathcal{N}[u_l]\,\rho_r\,dx}{\int_\Omega \bigl(\mathcal{N}[u_l]\bigr)^2 \rho_r\,dx + \gamma \int_{\partial\Omega} u_l^2\,\rho_b\,dx}\]

This is a least-squares projection: \(\alpha\) is chosen so that \(u_l\) corrects whatever systematic low-frequency residual \(u_h\) leaves behind.

8.7 Experimental Results 📊

The gradient scaling analysis clarifies the mechanism behind high-frequency amplification. For Fourier modes in the loss, gradient magnitudes scale as:

\[G_m(k, c) \propto k^{2m} |c|\]

where \(m = 0\) for \(L^2\) regression and \(m = 2\) for PINN residual losses (\(\mathcal{N}[u] - f\) involves second derivatives). The \(k^4\) factor in PINN losses naturally amplifies high-frequency gradients — the cross-attention architecture is particularly compatible with this regime because the frequency tokens can be dynamically upweighted to exploit the amplification.

8.8 PyTorch Implementation 🛠️

The architecture has three separable components: (1) the multiscale token bank \(H(x)\), (2) the cross-attention block with additive frequency masking, and (3) the network that stacks blocks and reads out. We build each in turn, then run a toy experiment.

Multiscale Token Bank

import torch
import torch.nn as nn
import torch.nn.functional as F

class MultiscaleTokenBank(nn.Module):
    """H(x) ∈ R^{N_tok × d_q}: dyadic-scaled random Fourier feature bank."""

    def __init__(self, d_in: int, M_base: int, K: int, d_q: int, sigma: float = 1.0):
        super().__init__()
        M = M_base * (K + 1)
        assert M % d_q == 0, "M_base*(K+1) must be divisible by d_q"

        # Base frequencies sampled once; fixed (not trained).
        omega_base = torch.randn(M_base, d_in) / sigma     # (M_base, d_in)

        # Dyadic dilations: omega_{m,k} = 2^k * omega_m.
        # scales shape: (K+1, 1, 1) to broadcast over (M_base, d_in).
        scales = 2.0 ** torch.arange(K + 1)                # (K+1,)
        omega = (scales[:, None, None] * omega_base[None]   # (K+1, M_base, d_in)
                 ).reshape(M, d_in)                         # (M, d_in)
        self.register_buffer('omega', omega)

        self.phase    = nn.Parameter(torch.zeros(M))        # learnable phase offsets
        self.log_beta = nn.Parameter(torch.tensor(0.0))     # amplitude envelope scale

        self.M, self.d_q = M, d_q

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x : (B, d_in)  →  H : (B, N_tok, d_q)
        beta = self.log_beta.exp()
        amp  = torch.exp(-beta * self.omega.norm(dim=-1))   # (M,)  amplitude envelope
        proj = x @ self.omega.T + self.phase                # (B, M) raw projections
        feats = amp * torch.cos(proj) / self.M ** 0.5      # (B, M) normalised features
        return feats.view(x.shape[0], -1, self.d_q)         # (B, N_tok, d_q)

Cross-Attention Block

class CrossAttentionBlock(nn.Module):
    """
    One layer: query q ∈ R^{d_q} attends over H ∈ R^{N_tok × d_q}.
    Residual structure: q ← q + CA(q, H) + FF(q + CA(q, H)).
    """

    def __init__(self, d_q: int, n_toks: int):
        super().__init__()
        self.W_Q = nn.Linear(d_q, d_q, bias=False)
        self.W_K = nn.Linear(d_q, d_q, bias=False)
        self.W_V = nn.Linear(d_q, d_q, bias=False)
        self.ff  = nn.Linear(d_q, d_q)
        self.d_q = d_q

        # Additive logit mask η ∈ R^{n_toks}, one scalar per token position.
        # Initialise all to 0 (unmasked). Posterior tokens can be clamped ≤ 0
        # externally to enforce the coarse-to-fine curriculum.
        self.mask = nn.Parameter(torch.zeros(n_toks))

    def forward(self, q: torch.Tensor, H: torch.Tensor) -> torch.Tensor:
        # q : (B, d_q)       H : (B, N_tok, d_q)
        Q = self.W_Q(q).unsqueeze(1)                        # (B, 1,     d_q)
        K = self.W_K(H)                                     # (B, N_tok, d_q)
        V = self.W_V(H)                                     # (B, N_tok, d_q)

        logits = (Q @ K.transpose(-2, -1)) / self.d_q**0.5 # (B, 1, N_tok)
        logits = logits + self.mask                         # additive mask (broadcasts)
        attn   = logits.softmax(dim=-1)                     # (B, 1, N_tok)

        ca_out = (attn @ V).squeeze(1)                      # (B, d_q)
        q_tilde = q + ca_out                                # CA residual
        return q_tilde + F.relu(self.ff(q_tilde))           # FF residual

Full Network

class RFFCA(nn.Module):
    """
    RFF-CA: multiscale token bank + L cross-attention blocks + linear readout.
    Default hyperparameters follow Feng et al. (2025) scaled down for a 1D toy.
    """

    def __init__(
        self,
        d_in:   int   = 1,
        M_base: int   = 32,   # base frequencies per octave
        K:      int   = 3,    # octaves (0..K); total M = M_base*(K+1)
        d_q:    int   = 32,   # token / query width
        L:      int   = 4,    # number of cross-attention blocks
        sigma:  float = 1.0,  # base frequency scale
        d_out:  int   = 1,
    ):
        super().__init__()
        N_tok = M_base * (K + 1) // d_q
        self.bank   = MultiscaleTokenBank(d_in, M_base, K, d_q, sigma)
        self.q0     = nn.Linear(d_in, d_q)             # initial query from raw coords
        self.blocks = nn.ModuleList([
            CrossAttentionBlock(d_q, N_tok) for _ in range(L)
        ])
        self.out = nn.Linear(d_q, d_out)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        H = self.bank(x)                    # (B, N_tok, d_q)
        q = self.q0(x)                      # (B, d_q) — query from coordinates
        for block in self.blocks:
            q = block(q, H)                 # each block refines q by attending H
        return self.out(q)                  # (B, d_out)

Toy Experiment: 1D High-Frequency Regression

Target: \(f^*(x) = \sin(2\pi \cdot 20 \cdot x)\) on \([0, 1]\) — frequency 20, well outside the effective range of a plain MLP (recall \(\omega^*_\text{plain} \approx T^{1/(d+1)} = T\) for \(d=1\); even \(T=10^4\) steps gives \(\omega^* \approx 10^4\), but with realistic width and learning rate the actual range is much lower).

import torch, torch.optim as optim

# ── data ───────────────────────────────────────────────────────────────────
torch.manual_seed(0)
N = 512
x_train = torch.linspace(0, 1, N).unsqueeze(-1)       # (N, 1)
y_train = torch.sin(2 * torch.pi * 20 * x_train)      # (N, 1)

# ── baselines ──────────────────────────────────────────────────────────────
class PlainMLP(nn.Module):
    def __init__(self, width=128, depth=4):
        super().__init__()
        layers = [nn.Linear(1, width), nn.ReLU()]
        for _ in range(depth - 2):
            layers += [nn.Linear(width, width), nn.ReLU()]
        layers.append(nn.Linear(width, 1))
        self.net = nn.Sequential(*layers)
    def forward(self, x): return self.net(x)

class RFFNN(nn.Module):
    """Fixed random Fourier features, no cross-attention."""
    def __init__(self, M=128, sigma=20.0, width=128, depth=3):
        super().__init__()
        omega = torch.randn(M, 1) * sigma
        self.register_buffer('omega', omega)
        self.phase = nn.Parameter(torch.zeros(M))
        layers = [nn.Linear(M, width), nn.ReLU()]
        for _ in range(depth - 2):
            layers += [nn.Linear(width, width), nn.ReLU()]
        layers.append(nn.Linear(width, 1))
        self.net = nn.Sequential(*layers)
    def forward(self, x):
        feats = torch.cos(x @ self.omega.T + self.phase) / self.omega.shape[0]**0.5
        return self.net(feats)

# ── models ─────────────────────────────────────────────────────────────────
models = {
    "PlainMLP": PlainMLP(),
    "RFF-NN":   RFFNN(sigma=20.0),               # hand-tuned sigma for fair comparison
    "RFF-CA":   RFFCA(d_in=1, M_base=32, K=3,   # sigma=5: octaves 5,10,20,40
                      d_q=32, L=4, sigma=5.0),
}

# ── training loop ──────────────────────────────────────────────────────────
STEPS = 5_000
losses = {name: [] for name in models}

for name, model in models.items():
    opt = optim.Adam(model.parameters(), lr=1e-3)
    for step in range(STEPS):
        pred = model(x_train)
        loss = F.mse_loss(pred, y_train)
        opt.zero_grad(); loss.backward(); opt.step()
        losses[name].append(loss.item())
    print(f"{name:10s}  final L2 = {losses[name][-1]:.4e}")

Expected output (approximate — varies by seed):

PlainMLP    final L2 = 2.3e-01   # stuck at ~half-variance: spectral bias
RFF-NN      final L2 = 8.1e-03   # better, but sensitive to sigma choice
RFF-CA      final L2 = 1.4e-05   # near-perfect: attention routes to freq-20 tokens

Inspecting the Attention Weights

After training, the learned attention map reveals which frequency tokens each block routes to:

model = models["RFF-CA"]
model.eval()
with torch.no_grad():
    x_probe = torch.tensor([0.25](/notes/0.25/))       # a single probe point
    H = model.bank(x_probe)               # (1, N_tok, d_q)
    q = model.q0(x_probe)

    for i, block in enumerate(model.blocks):
        Q = block.W_Q(q).unsqueeze(1)
        K = block.W_K(H)
        logits = (Q @ K.transpose(-2, -1)) / block.d_q**0.5 + block.mask
        attn = logits.softmax(dim=-1).squeeze()   # (N_tok,)
        top = attn.argmax().item()
        print(f"Block {i}: peak attention on token {top}  "
              f"(weight {attn[top]:.3f})")
        q = block(q, H)

In a trained model, later blocks show increasingly concentrated attention — early blocks spread across octaves while the final blocks collapse onto the frequency-20 tokens. This is the learned coarse-to-fine curriculum emerging from gradient descent rather than the hand-scheduled masking of §8.4.

9. References