Conceptio › Archive › arXiv CS

arXiv CSopen access

Learning Orthonormal Bases for Function Spaces

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

neuralnetworks

machine learning, deep learning, neural networks

Learning Orthonormal Bases for Function Spaces

arXiv:2605.19959v1 [cs.LG] 19 May 2026

Hamidreza Kamkari MIT CSAIL Cambridge, MA [email protected]

Mohammad Sina Nabizadeh MIT CSAIL Cambridge, MA [email protected] Justin Solomon MIT CSAIL Cambridge, MA [email protected]

Abstract Infinite-dimensional orthonormal basis expansions play a central role in representing and computing with function spaces due to their favorable linear algebraic properties. However, common bases such as Fourier or wavelets are fixed and do not adapt to the structure of a given problem or dataset. In this paper, we aim to represent these bases with neural networks and optimize them. Our key idea is that any target infinite-dimensional orthonormal basis can be viewed either as a point on the Lie manifold of the orthogonal group, or equivalently, as the endpoint of a continuous path on that manifold that connects a reference basis, e.g. Fourier, to that target. Paths on the Lie manifold satisfy ordinary differential equations (ODEs) governed by skew-adjoint integral operators. Using neural networks to define finite-rank generators of such ODEs allows us to parameterize and optimize orthonormal bases in function space. While relying on finite-rank generators to model infinite operators might seem restrictive, we prove a universality result: even with a rank-2 generator, the integrated solutions of the ODE are dense in the orthogonal group under the appropriate operator topology. In other words, for any target orthonormal basis, there exists a path originating from a reference basis and driven by finite-rank generators that gets arbitrarily close to that target basis. We demonstrate the flexibility of our framework by transforming the Fourier basis into the principal components of a functional dataset, eigenfunctions of linear operators, or dynamic modes of energy-preserving physical simulations.

Introduction

Orthonormal basis (ONB) expansions are basic building blocks of various fields. They help us represent functions, run numerical simulations, and analyze data in machine learning. The choice of basis governs not only approximation accuracy but also the computational cost of downstream tasks. Often, the most useful basis is adapted to the problem at hand. For instance, in data analysis, principal components optimally capture variance, allowing one to analyze the data by looking at only the first few basis elements [50, 1]. In physics and geometry, the first few spectral eigenfunctions approximate the dynamic or geometric modes of the system [10, 46]. In signal processing, Fourier and wavelet series [39, 6] help pinpoint various signal characteristics. Neural networks flexibly parameterize complex function spaces and adapt to specific tasks via gradient-based optimization. Here, we leverage their expressiveness for the space of ONBs and ask: Can we parameterize infinite series of orthonormal and continuous bases with neural networks? Preprint.

Transformed basis element

principal component analysis

optimizing

diagonalizing linear operators

4 Fourier basis

dynamic modes of ﬂuids (a) Basis adaptation across tasks

(b) Discretization-free property

Figure 1: Method and Applications Overview. (a) Parameterizing a change-of-basis map Qθ allows us to adapt a reference ONB to various applications: eigenfaces that describe the variance of the data, visualizing the training dynamics of a non-linear neural function by diagonalizing its neural tangent kernel, or dynamic modes of fluid-flow. (b) The learned bases inherit properties of the reference such as being discretization-free. (c) We model paths on the Lie manifold of the infinite orthogonal group and gradually transform an initial basis by traversing this manifold; by optimizing the path, we are able to learn ONBs.

Classical ONBs such as Fourier series provide a good starting point: they are infinite and orthonormal by construction, inner products decompose cleanly into coefficient sums, and projection requires a single inner product. They are also continuous and discretization-free: they represent functions via linear combination of other functions on the entire domain, rather than as discrete values on a grid. Thus, the resulting space can be evaluated, differentiated, or integrated at any point in the domain. Any two ONBs, up to sign, are connected by a special orthogonal change-of-basis operator. Therefore, rather than directly modeling the space of ONBs, we can parameterize the equivalent orthogonal map that inputs a reference ONB—e.g. the Fourier series—and outputs a basis tailored to the application at hand. This construction allows the target ONB to inherit desirable properties of this reference: remaining infinite, continuous, and orthonormal. Figure 1 (a) shows how parameterizing this operator allows us to map the Fourier basis to the principal components of the celebrity face dataset [29], recover the spectrum of a linear operator, or extract the dynamic modes of a system. To tractably parameterize this change-of-basis operator, we draw on Lie theory [24]. Rather than optimizing an orthogonal operator directly, we view it as a continuous path across the manifold of the orthogonal group (see Figure 1; c). Geometrically, this path corresponds to a sequence of infinitesimal rotations in function space, modeled as the flow of an ODE governed by a timedependent, skew-adjoint infinitesimal generator. The orthogonal group forms a highly constrained manifold that is hard to parameterize, but the space of skew-adjoint generators is linear and easier to model with a neural network. The remaining challenge is that even the space of generators is infinite-dimensional. To resolve this, we prove in section 4 that any path on this infinite group can be arbitrarily well-approximated using a finite-rank generator. Ultimately, this reduces the task of parameterizing and optimizing an infinite ONB to the following: we use a neural network to model a finite-rank, skew-adjoint generator and integrate the resulting ODE to smoothly change a reference basis into a target basis. By optimizing the network’s parameters via neural ODE backpropagation [9], we model a continuous path on the orthogonal group that adapts to our needs. Our contributions are: (i) reformulating orthonormal bases as a rotation path on the infinitedimensional orthogonal group, achieved by solving an ODE; (ii) a universality theorem showing that a large family of rotations (specifically, those generated by Hilbert-Schmidt operators) can be approximated by the flow of an ODE governed by a low-rank skew-adjoint operator; (iii) a practical neural parameterization of this generator, trained end-to-end via ODE-based backpropagation, yielding an orthonormal, continuous, and infinite basis that can be adapted to variational problems; (iv) a principled variational objective for diagonalizing linear operators, demonstrated in two applications by recovering the eigenfunctions of the data covariance operator and the neural tangent kernel; and (v) applications beyond eigenproblems to modeling dynamical systems, demonstrating that our basis transformations can approximate the Koopman operator of fluids.

Related Work

Basis expansions. A principled way to discover a problem-adapted basis is to compute the eigenfunctions of an underlying linear operator. In physics, eigenfunctions of differential operators—such 2

as spherical harmonics, Bessel functions, and Hermite polynomials—diagonalize the Laplacian and form the foundation of spectral methods for solving differential equations [23, 53, 34, 8]. Alternatively, operators can be data-driven rather than analytic; for example, in principal component analysis (PCA), the principal components of a stochastic process are the eigenfunctions of its covariance operator [19, 1]. In other cases, bases are engineered rather than linked to operators; wavelet bases are one such example which are designed to capture local signal characteristics [39]. Despite their different origins, ONBs are useful and unanimous across various fields involving numerical processing. Eigendecomposition. Solving an eigenproblem is a standard method for obtaining custom bases. The classical approach involves discretizing the operator as a matrix on a grid or mesh and applying numerical eigensolvers [45, 22, 50]. However, these approaches inevitably incur discretization error in continuous state spaces. While solving the discretized system and applying post-hoc interpolation is a possible workaround [5], recent work instead suggests using neural networks to learn continuous eigenfunctions directly by solving variational problems. For instance, Spectral Inference Networks [43] pair this approach with an explicit orthogonalization step, whereas NeuralEF [12] avoids it by extending EigenGame [18] to function space. Further extending this paradigm, Ben-Shaul et al. [4] target differential self-adjoint operators on free-form domains, and Azmoodeh et al. [3] address continuous-time PCA via the directly representing the ONBs. While these neural methods are discretization-free, they typically recover only a small, finite number of eigenfunctions without strictly enforcing orthogonality. Our framework, by contrast, maintains exact orthonormality across the full eigenfunction sequence by always remaining on the Lie manifold. Non-linear representations of function spaces. An alternative to using orthonormal bases is to forego the linear structure of function spaces entirely in favor of highly expressive, non-linear neural parameterizations. Neural operators [35, 36, 38, 31, 37], implicit neural representations (INRs) [48, 40, 55, 14], and physics-informed neural networks (PINNs) [44] loosely fall into this category and have achieved success across image processing, computer graphics, and physical simulation [41, 47]. However, while these representations adapt remarkably well to specific problems, they sacrifice the linear algebraic structure of the function space: linear operations in function space— such as addition or applying a differential operator—map to highly non-linear and unpredictable transformations in the network’s weight space [56, 57, 11]. While recent architectures like the Fourier tensor network [2] have begun to address this by, for example, explicitly structuring INRs with orthogonal Fourier bases and low-rank tensor coefficients, our framework offers an alternative: we leverage the adaptive power of neural networks while strictly preserving the linear algebraic structure of the function space by parameterizing bases. We even demonstrate how to bridge these two worlds; in one of our applications, we take an arbitrary non-linear neural network and project it back into a structured function space by learning the continuous ONB that diagonalizes its neural tangent kernel.

Mathematical Preliminaries

Function spaces. We work with functions defined on a d-dimensional compact domain Ω equipped 2 with finite measure µ. We consider the real-valued Hilbert R space H := L (Ω, µ), with its inner product and norm defined in standard fashion: ⟨f, g⟩ := Ω f (ω)g(ω) µ(dω) = Eω∼µ [f (ω)g(ω)] and ∥f ∥2 := ⟨f, f ⟩. Since µ is finite, we may without loss of generality assume it is a probability measure, allowing for Monte Carlo estimation of the inner product. We assume H is separable, meaning it admits a countable complete ONB {φn }∞ n=1 . Concretely, for any f ∈ H, limn→∞ ∥f − Pn ⟨f, φm ⟩ · φm ∥ = 0. For simplicity, we assume Ω possesses a standard Euclidean metric; thus, m=1 a valid ONB {φn }∞ n=1 (among many) can be the standard d-dimensional real Fourier basis. This basis may be suboptimal for a given problem, motivating us to change it to a problem-adapted one. Henceforth, we will use {φi }∞ i=1 to denote the initial (reference) basis. Linear operators. A bounded linear transformation on H is a linear operator. Any such operator A admits a unique adjoint denoted by A∗ , defined such that: ⟨f, Ag⟩ = ⟨A∗ f, g⟩, ∀f, g ∈ H. When A can be written in an integral kernel form A[·, ·] : Ω × Ω → R, its action follows the identity: Z (Af )(x) = A[x, y] f (y) µ(dy); consequently, ∀x, y ∈ Ω : A∗ [x, y] = A[y, x]. (1) Ω

A is skew-adjoint if A∗ = −A, or equivalently, ∀x, y ∈ Ω : A[x, y] = −A[y, x]. 3

The space of operators themselves have norms and topology. Here, we use the Hilbert-Schmidt norm: ∥A∥2HS :=

∞ X

∥Aφn ∥ =

n=1

∞ X

⟨Aφn , Aφn ⟩ =

n=1

∞ X

⟨φn , A∗ Aφn ⟩ =: tr(A∗ A).

(2)

n=1

where {φn }∞ n=1 can be any complete ONB. An operator A is Hilbert-Schmidt if this norm is finite—or equivalently, if A∗ A is trace-class, meaning its trace is finite and well-defined. ∞ Infinite-dimensional orthogonal group. Any two ONBs, {φi }∞ i=1 and {ϕi }i=1 , are related by an ∗ orthogonal change-of-basis operator Q : H → H, where Q Q = I.

To make our search for a basis tractable, we borrow a concepts from Lie theory [24]: a large subset of orthogonal operators Q can be characterized by a smooth manifold obtained via the exponentiation of a skew-adjoint operator. Concretely, let Q(t) be a family of operators obeying the following ODE: d Q(t) = KQ(t), dt

where K : H → H is skew-adjoint.

Q(0) = I,

(3)

The solution, Q(t), remains orthogonal because the inner product Q∗ Q is preserved: d ∗ (4) (Q Q) = Q̇∗ Q + Q∗ Q̇ = Q∗ K ∗ Q + Q∗ KQ = Q∗ (K ∗ + K)Q = 0. dt Searching over the smooth family of generators K instead of Q, establishes a tractable search space. Formally, this allows us to work with the following continuous group: Definition 1. The Hilbert-Schmidt orthogonal group, denoted SO(HS), is the space of operators that can be written as solutions to the ODE in (3), where K is additionally Hilbert-Schmidt.

Approximation Theory of the Orthogonal Group

Our goal is to represent an orthogonal change-of-basis operator Q ∈ SO(HS) using the parameters θ of a neural network. As a first step, rather than directly modeling Q, we model its skew-adjoint infinitesimal generator Kθ (t) as a time-evolving operator: d Q(t) = Kθ (t)Q(t), dt

Q(0) = I,

(5)

integrated to time T > 0. Even when Kθ (t) is time-dependent, (4) shows Q(t) is an orthogonal time-evolving operator. In a slight abuse of notation, we will denote this operator by Qθ := Q(T ), dropping the time evaluation and using θ emphasize its parameterization. We have reduced the problem of designing an orthogonal operator to designing a skew-adjoint operator with a neural network. If the infinitesimal generator Kθ were to have finite rank r = 2, we could represent its skew-adjoint kernel using two unrestricted multilayer perceptrons (MLPs): Kθ (t)[x, y] = MLP1 (x, t; θ) · MLP2 (y, t; θ) − MLP2 (x, t; θ) · MLP1 (y, t; θ).

(6)

Although designing the generator to have a finite rank might appear as a restrictive choice, integrating low-rank dynamics over time allows one to approximate any transformation in SO(HS). This is true even when the rank is as small as r = 2. We now make this concrete: Theorem 1 (Approximating the Orthogonal Group via Rank-2 Generators). Let QT denote the set of all operators Q(T ) obtained by solving the initial value problem (5) at time T > 0, where the generator K(t) is uniformly Hilbert-Schmidt, continuous, and rank-2. Then, QT is dense in SO(HS) with respect to the Hilbert-Schmidt norm topology. In other words, every rotation in this infinite orthogonal group is arbitrarily well-approximated by some time-evolving rank-2 generator. Proof Sketch. We first approximate an arbitrary Q ∈ SO(HS) with an operator that acts non-trivially only on a finite n-dimensional subspace. Restricted to this subspace, we have a constructive proof: this operator is isomorphic to an n × n orthogonal matrix. A classical result from linear algebra states that any finite-dimensional rotation can be factorized into a sequence of planar Givens’ rotations [22]. Each planar rotation is generated by a fixed rank-2 skew-adjoint operator and by using a timedependent scalar control function, we can chain these stationary generators in sequence; effectively 4

“turning on” each planar rotation one after the other. This constructs a single continuous time-varying rank-2 generator whose time-ordered integral results in the target finite-dimensional rotation. Finally, taking the limit of n → ∞ in the Hilbert-Schmidt operator topology extends this approximation to the entire infinite-dimensional group SO(HS); for full details, see Appendix A. Theorem 1 thus establishes that an entire change-of-basis can be implicitly parameterized through a finite-rank Kθ as in (6). Due to the universal approximation theorem of MLPs (the version that uses compactness of Ω by Hornik et al. [26]), even the rank r = 2 model in (6) with a sufficiently overparameterized θ can implicitly represent an arbitrary basis in function space.

Variational Problems: Adapting the Basis to Various Tasks

The previous sections establish how a path on the orthogonal group manifold can be parameterized with a neural network and traversed while preserving structure. Our next task is to put this parameterization to use. To that end, we derive variational objective functions that reward or penalize paths on the Lie manifold and allow us to optimize it and find useful bases. Eigenproblems. One class of variational problems is finding Qθ that diagonalizes an infinitedimensional linear operators A : H → H; this is akin to finding eigenvectors of n × n matrices in the finite case. In our formulation, we can freely search the space of operators Q ∈ SO(HS) to find one that diagonalizes A. To do so, we reformulate the eigenproblem in the variational language: ∞ Theorem 2. Let A be a trace-class positive semidefinite linear P∞ operator, and let {p(i)}i=1 be a strictly decreasing sequence of positive weights such that i=1 p(i) = 1, and {φi }∞ i=1 be any orthonormal basis on H. Consider Qopt as the solution to the following variational problem: X max p(i) · ⟨Qφi , AQφi ⟩ (7) Q∈SO(HS) i=1

Then, ϕ⋆i := Qopt φi are the eigenfunctions of A. That is, Aϕ⋆i λ1 ≥ λ2 ≥ λ3 ≥ · · · are the corresponding eigenvalues.

= λi ϕ⋆i for all i ≥ 1, where

See Appendix B for a proof. Interpreting p(·) as a probability distribution on indices i allows us to write the objective as an expectation. Furthermore, replacing Q with Qθ changes an otherwise intractable search over SO(HS) into one over parameters θ, compatible with stochastic optimization: max θ

Ei∼p ⟨Qθ φi , AQθ φi ⟩.

(8)

Substituting various operators for A allows us to solve different eigenproblems in function space. We provide an example below; Appendix B details another involving the neural tangent kernel. PCA in function space. We formulate PCA in function space as an eigenproblem.1 In this setting, we assume our data consists of random functions X drawn from an Ω-indexed stochastic process. This process can be formally viewed as a probability distribution P over the Hilbert space H. For simplicity of exposition, we assume that these functions have zero mean, EX∼P [X] = 0; we defer the treatment of the non-centered case to Appendix B. The covariance operator C : H → H is defined by its action on any function f ∈ H as Cf = EX∼P [X⟨X, f ⟩]. Plugging this definition into (8) yields the following derivation: (9) Ei∼p ⟨Qθ φi , CQθ φi ⟩ = Ei∼p Qθ φi , EX∼P X⟨X, Qθ φi ⟩ = Ei∼p,X∼P ⟨X, Qθ φi ⟩2 . This objective maximizes the explained variance of the stochastic process P. Practically, the optimization algorithm is as follows: sample i ∼ p, retrieve the corresponding initial (Fourier) basis function φi , apply Qθ , sample a random function X ∼ P, estimate the squared projection ⟨X, Qθ φi ⟩2 using b finally maximize via stochastic gradient ascent on θ. random quadrature points Ω, Koopman operators. We also cover an application that cannot be formulated as an eigenproblem. For a dynamical system governed by a state transformation Ψ : Ω → Ω, the associated Koopman operator U [30] acts on observable functions via composition Uf := f ◦ Ψ. When Ψ is volume-preserving (e.g. incompressible flow), U becomes orthogonal. Because Qθ is orthogonal by design, it provides a 1 PCA in function space is sometimes referred to as the Karhunen-Loéve expansion in the literature [49, 20, 54].

Solver exploded Forward Euler Backward Euler Cayley Transform

Basis

Transformed Fourier basis

energy

Neural network outputs

Initial Fourier basis

Random quadrature points Random time discretization (a) Structure-preserving Cayley Integration

Time

(b) Ablating various solvers

Figure 2: Numerical Pipeline. (a) A Fourier basis is gradually transformed into principal functions of CelebA dataset. We b = {ωj }D ∼ µ and time-discretization {t0 = 0, . . . , tL = T } and perform L Cayley steps. select quadrature points Ω j=1 Each step takes O(r3 + r2 D) that exploits the rank r structure of the generator. (b) Using off-the-shelf ODE solvers result in loss or explosion of the function norm, by contrast, Cayley integration keeps it intact.

suitable inductive bias to learn these dynamics. To that end, we directly fit Qθ to the Koopman U by minimizing the expected projection error: min θ

Ei∼p ∥Qθ φi − φi ◦ Ψ∥2 .

(10)

Since both Qθ and U are linear, training them to match on the complete basis {φi } is sufficient to fit Qθ to U. In turn, the learned ONB ϕi = Qθ φi represents the dynamic modes of the system. This formulation becomes particularly valuable when Ψ cannot be numerically simulated in a strictly volume-preserving manner. In such cases, learning Qθ that is orthogonal by design guarantees energy preservation in the system—∀f ∈ H : ∥Qθ f ∥2 = ∥f ∥2 ; even when ∥f ◦ Ψ∥2 ̸= ∥f ∥2 numerically. In Figure 6, we demonstrate this via a simple fluid dynamical flow: by repeatedly transforming a Fourier basis using Qθ , we generate future state snapshots that strictly preserve norm, despite being trained on imperfect numerical approximations of Ψ.

Algorithmic Details

Our objectives require the action of Qθ on the initial basis functions φi . Defining ϕθi (t) := Qθ (t)φi and substituting into (5), we obtain d θ ϕ (t) = Kθ (t)ϕθi (t), ϕθi (0) = φi . (11) dt i In principle, one could integrate (11) with off-the-shelf ODE solvers and optimize θ. However, standard solvers like forward and backward Euler accumulate numerical error over time and violate the orthonormal structure; we have included a comparison in Figure 2 with more details in Appendix C. We therefore use a structure-preserving Lie group integrator [13, 27, 7] that preserves inner products. Structure-preserving ODE integration. Our Lie integrator is the Cayley method [27], a structure preserving method which discretizes time as t0 = 0 ≤ t1 < . . . < tL = T and advance the system ∆tℓ b from tℓ−1 to tℓ . For each ℓ, the generator is approximated by a midpoint rule Kℓ = Kθ (tℓ−1 + 2 ) b ℓ . An important step in the Cayley method and the system is evolved by applying Qℓ = exp ∆t · K is that it replaces the exponential with its (1, 1)-Padé approximant which preserves inner products: −1 b ℓ ϕθi (tℓ−1 ), b ℓ := I − ∆t K bℓ b ϕθi (tℓ ) ≈ Q Q I + ∆t (12) 2 2 Kℓ ≈ Qℓ . b ℓ is skew-adjoint, then Q b ℓ is orthogonal [22, 17], so the Cayley One can verify algebraically that if K update preserves the norms and inner products of the basis exactly at each step. b = {ωj }D ∼ µ, where the ωj Applying (12) requires discretization at some quadrature points Ω j=1 are drawn via stratified sampling. To avoid overfitting to any discretization, we resample the points b ℓ ∈ RD×D , and K b ℓ ∈ RD×D are no before each ODE solve. Hence, discretely, ϕbθi (t) ∈ RD , Q longer functions and operators, but random vectors and matrices. Similarly, we resample in each training iteration tℓ ∼ Unif(0, T ) to avoid overfitting to any particular time discretization. Low-rank parameterization. While it is mathematically valid to parameterize Kθ (t) using a rank-2 generator in (6), we generalize it to higher ranks r ≥ 2 for more expressiveness. In particular, Kθ (t) = Uθ∗ (t) Mθ (t) − Mθ (t)⊤ Uθ (t), (13) 6

Fourier

PCA

Ours

(a) Samples of 1D signals

(b) Comparison of Bases

max cutoﬀ freq. of 16 max cutoﬀ freq. of 32

Ours

Fourier

PCA

Ours

Reconstruction error

max cutoﬀ freq. of 4

Fourier

max cutoﬀ freq. of 2

Fourier basis Learned basis (Ours) PCA number of basis used

PCA

(d) Inset comparisons

(e) Residual energy of reconstruction

Figure 3: 1D Synthetic PCA on Function Space. (a) Samples from a distribution of signals described by (15), exhibiting a jump at x = 0.5 on a 1D domain. (b) The Fourier basis is transformed to diagonalize the covariance operator, allowing the basis to better capture the variance of signals. Our results improve upon the traditional finite-dimensional PCA construction, which is limited by the Nyquist cutoff inherent to that setting. (c) Reconstruction of several sample signals from the distribution. Our learned basis yields better reconstructions across all cutoff frequencies. (d) In particular, finite PCA is limited by Nyquist limits and the Fourier basis exhibits the Gibbs ringing phenomenon at the discontinuity. (e) A quantitative comparison of reconstructions averaged across the distribution of signals. The reconstruction error is shown on a log scale. We achieve an order-of-magnitude improvement using roughly 100 basis elements.

where Uθ (t) : H → Rr has r rows (i.e. r functions in H) and Mθ (t) ∈ Rr×r is an arbitrary matrix. By construction, Kθ (t) is skew-adjoint, has rank ≤ r, and depends smoothly on θ. In practice, Uθ is implemented as a network that inputs both time t and coordinate ω ∈ Ω and outputs r values, and Mθ is implemented as another network that inputs time t and outputs an r × r matrix. Numerically, ℓ b b b b∗ Kθ is only evaluated on the midpoints tℓ + ∆t 2 and quadrature points Ω, producing Kℓ , Uℓ , and Uℓ . Efficient Solver. The quadrature set must be large enough (1024 ≤ D ≤ 4096 for most experiments) b ℓ to ϕbθ to accurately approximate the underlying function-space operations. Naïvely applying Q i 3 would require O(D ) matrix vector products and linear solves; however, leveraging the low-rank b ℓ allows us to be much more efficient: Each step amounts to applying Q b ℓ to the input ϕbθ . structure of K i b ℓ can be done in O(r2 + rD) Due to the low-rank structure, the forward action of applying I + 12 K b ℓ )−1 , we apply the Woodbury identity [25]: time. On the other hand, for the inverse action (I − 12 K

bℓ I − 21 K

−1

−1 bℓ∗ Sbℓ Ir×r − 1 U bℓ U bℓ∗ Sbℓ bℓ , = I + 12 U U 2

(14)

∆tℓ ⊤ ∗ r×r ℓ where Sbℓ := Mθ (tℓ + ∆t is a small matrix of 2 ) − Mθ (tℓ + 2 ) . The Gram matrix Uθ Uθ ∈ R b with O(r2 D) operations. Hence, the total backward action inner products, which we estimate on Ω reduces to O(r3 + r2 D) linear solves and matrix vector products, which are all linear in D.

Figure 2 illustrates the algorithm in action. For full details on the integrator, the values of L, r, and D, ablations, and the architecture used for each component, see Appendix C. The algorithm presented here only involves real-valued functions, for experiments involving multiple channels, such as the CelebA dataset with RGB, we use a multichannel variant explained in the same appendix.

Experiments

In this section, we present experiments that evaluate our learned bases. Each task requires specific configuration choices. We focus on the primary takeaways and defer details and additional experiments to Appendix D. Our codebase is available at [https://github.com/HamidrezaKmK/Learning-ONBs]. 7

(b) CelebA Image dataset Cumulative captured energy

Learned basis (Ours) Original

Learned Fourier basis (Ours) basis

Fourier basis Original

(a) MNIST INR (SIREN zoo) dataset

Learned (Ours)

number of basis elements number of basis elements

(d) MNIST INR Cumulative captured energy

number of basis elements

(e) CelebA Cumulative captured energy

Figure 4: 2D PCA in Function Space. (a) Fourier basis transformed to explain the INR zoo of SIREN networks representing the MNIST dataset. Note that the resulting basis represents a dataset that is itself discretization-free. (b) Transformed bases for the CelebA dataset. (c) Reconstruction of several datapoints using the Fourier basis and our transformed, optimized basis. With only a few basis elements, our basis reconstructs the signals more faithfully and captures the shape of the digits better than the Fourier baseline. (d) Quantitative comparison of the cumulative captured energy (mean ± std across the dataset). Our learned basis captures more energy than the Fourier baseline. Because our basis is infinite, it gradually achieves full reconstruction without a Nyquist limit; any residual error stems from numerical integration.

1D functional PCA. We evaluate our approach using a synthetic 1D experiment on the domain Ω = [0, 1] equipped with a uniform measure µ. To create a deliberately challenging scenario for standard Fourier methods, we sample functions with a sharp discontinuity at x = 0.5 (see Figure 3): sin(2πkx) + 2x if x ≤ 0.5, f (x) = σ (15) − sin(2πkx) + 2x − 2 if x > 0.5, where the scale factor σ ∼ Rademacher takes values in {−1, 1} with equal probability, and the frequency multiplier k ∼ Geometric(1/3) is a strictly positive integer. While discontinuous signals are notoriously difficult for standard Fourier reconstruction, our method overcomes this issue, also surpassing the baseline finite PCA reconstruction with as few as 10 basis functions. Moreover, as shown in Figure 3, while finite PCA is optimal for discrete data, its grid-based discretization prevents it from fitting signals beyond the Nyquist rate. Our formulation bypasses this bottleneck. Real Data PCA. To maintain a discretization-free setting, we use the INR Zoo [42], an MNIST [33] dataset where each datapoint is a SIREN network [48] overfitted to an image. We have a similar setup for a subset of 1024 images from CelebA, using quadrature points on a 64 × 64 grid since INRs were unavailable. In both cases, our method learns to transform the Fourier basis into principal functions of these datasets in Figure 4. As illustrated, the learned basis achieves better reconstruction than the Fourier baseline and the first few basis elements capture most of the signal energy of the dataset. Diagonalizing the neural tangent kernel (NTK). Moving beyond PCA, we diagonalize the NTK operator at various stages of training a two-moon classifier. A short summary on what diagonalization of the NTK means is provided in Appendix B. Figure 5 compares against a grid-based eigensolver. While the discrete solver isolates eigenfunctions with lower eigenvalues slightly better, its fixed resolution fails to capture fine spatial details, particularly at the sharp center of the decision boundary. Our discretization-free method successfully resolves these localized features, highlighting a fundamental advantage of solving eigenproblems entirely in function space. Koopman operators. We train Qθ to approximate the Koopman operator, by optimizing the objective in (10). Here, the underlying dynamical system Ψ is the Taylor-Green vortex on a flat torus [52], simulated for a discrete time step of ∆t = 0.5s. Specifically, the discrete flow map Ψ is obtained by integrating the governing ODE ẋ = v(x) over ∆t using a fourth-order Runge-Kutta (RK4) scheme, where v parameterizes the Taylor-Green vortices. After training, we evaluate the model on a sample initial condition (see Figure 6) by iteratively applying the learned operator to generate extended fluid-flow trajectory snapshots. We compare the energy conservation of these predicted rollouts against various solvers. Although baselines that run Ψ may produce more visually plausible dynamics, they fail to preserve energy, evaluated as the Hilbert norm of the signal. In contrast, despite being noisier, our learned Koopman operator yields plausible results and strictly conserves energy. 8

inset

step 5000

step 500

step 250

step 1

inset

(a) NTK

(b) NTK diagonalized with grid solver

RK2

Ours

energy

Figure 5: Two-moon NTK. Diagonalizing the NTK over the course of training a two-moon classifier (steps 1, 250, 500, and 5000). (a) The evolving logit landscape and training data. (b) Top eigenfunctions recovered by a discrete grid-solver; the grid’s fixed resolution leads to noticeable pixelation and a loss of high-frequency information. (c) Eigenfunctions obtained using our proposed discretization-free method. Our approach matches the global structure of the leading eigenfunctions while resolving the fine, intricate details near the decision boundary. Note that the eigenspaces are comparable up to sign which explains why the colormaps sometimes appear swapped between the two methods.

RK4

Learned transport (Ours) RK1 transport RK2 transport RK4 transport Ground truth

Time (s)

(a) Advection in time of Taylor-Green Vortex

(b) Function energy over time

Figure 6: Learning the Koopman Operator. (a) Visual comparison of fluid-flow trajectories generated by our learned Koopman operator vs. the composition operator integrated via RK methods; our method is noisier but conserves energy. (b) L2 energy tracking over time, comparing standard solvers (1st , 2nd , and 4th -order RK) against iterative rollouts from Qθ .

Conclusions, Limitations, & Future Work

In this work, we introduced a neural framework for learning continuous orthonormal bases in function spaces. Our approach allows us to optimize these bases for various tasks while preserving their essential linear-algebraic properties, such as continuity and orthonormality. The Lie ODE formulation allowed us to prove theoretical guarantees on the expressiveness of these networks. However, optimizing and backpropagating through any ODE trajectory is well-known to be practically challenging [21, 16]. Reaching our theoretical limits often demands long integration horizons, which can complicate optimization dynamics and increase training times. This motivates future work on tailored architectures and optimizers to scale our method. Several exciting directions remain for future research. Here, we showed how one can diagonalize the covariance and NTK operators; a natural extension is to diagonalize elliptic operators in physics or the Laplace-Beltrami operator in differential geometry. Additionally, we showed how our change-of-basis operator can match the Koopman operator, which is also an integral component of the more general family of dynamic mode decomposition (DMD) [46] methods. It is natural to extend our method in the context of DMD with more complicated physical systems and their Koopman operators. In conclusion, we hope that parameterizing orthonormal bases as a continuous, optimizable neural network will offer a powerful new lens for modeling complex systems involving function spaces. 9

Acknowledgments and Disclosure of Funding The MIT Geometric Data Processing Group acknowledges the generous support of National Science Foundation grants IIS2335492 and OAC2403239, from the CSAIL Future of Data and FinTechAI programs, from the MIT–IBM Watson AI Laboratory, from the Wistron Corporation, from the MIT Generative AI Impact Consortium, from the Toyota–CSAIL Joint Research Center, and from Schmidt Sciences. The authors acknowledge further support from SideFX Software. The authors additionally declare that they have no competing interests. We thank Clément Jambon, Artem Lukoianov, Christopher Scarvelis, Ana Dodik, Alice Petrov, Dolores Garcia Marti, and Ahmed Mahmoud for providing valuable feedback on our work. Hamidreza Kamkari extends a special thank you to Yeganeh Gharedaghi for her instrumental contributions during the early stages of this project and her invaluable help with ideation.

References [1] Hervé Abdi and Lynne J Williams. Principal component analysis. Wiley interdisciplinary reviews: computational statistics, 2(4):433–459, 2010. [2] Pooya Ashtari, Pourya Behmandpoor, Nikos Deligiannis, and Aleksandra Pizurica. Futon: Fourier tensor network for implicit neural representations. arXiv preprint arXiv: 2602.13414, 2026. [3] Shayan K Azmoodeh, Krishna Subramani, and Paris Smaragdis. Continuous-time signal decomposition: An implicit neural generalization of pca and ica. In 2025 IEEE 35th International Workshop on Machine Learning for Signal Processing (MLSP), pages 1–6. IEEE, 2025. [4] Ido Ben-Shaul, Leah Bar, Dalia Fishelov, and Nir Sochen. Deep learning solution of the eigenvalue problem for differential operators. Neural Computation, 35(6):1100–1134, 2023. [5] Yoshua Bengio, Jean-françcois Paiement, Pascal Vincent, Olivier Delalleau, Nicolas Roux, and Marie Ouimet. Out-of-sample extensions for lle, isomap, mds, eigenmaps, and spectral clustering. Advances in neural information processing systems, 16, 2003. [6] Pierre Brémaud. Mathematical principles of signal processing: Fourier and wavelet analysis. Springer, 2002. [7] Elena Celledoni, Håkon Marthinsen, and Brynjulf Owren. An introduction to lie group integrators–basics, new developments and applications. Journal of Computational Physics, 257: 1040–1061, 2014. [8] Yue Chang, Otman Benchekroun, Maurizio M. Chiaramonte, Peter Yichen Chen, and Eitan Grinspun. Shape space spectra. ACM Trans. Graph., 44(4), July 2025. ISSN 0730-0301. doi: 10.1145/3731148. URL https://doi.org/10.1145/3731148. [9] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018. [10] Richard Courant and David Hilbert. Methods of mathematical physics, vol. i. Phys. Today, 7 (5):17–17, 1954. [11] Luca De Luigi, Adriano Cardace, Riccardo Spezialetti, Pierluigi Zama Ramirez, Samuele Salti, and Luigi Di Stefano. Deep learning on implicit neural representations of shapes. The Eleventh International Conference on Learning Representations, 2023. [12] Zhijie Deng, Jiaxin Shi, and Jun Zhu. Neuralef: Deconstructing kernels by deep neural networks. In International Conference on Machine Learning, pages 4976–4992. PMLR, 2022. [13] Fasma Diele, Luciano Lopez, and R Peluso. The cayley transform in the numerical solution of unitary differential systems. Advances in computational mathematics, 8(4):317–334, 1998. 10

[14] Emilien Dupont, Hyunjik Kim, SM Eslami, Danilo Rezende, and Dan Rosenbaum. From data to functa: Your data point is a function and you can treat it like one. arXiv preprint arXiv:2201.12204, 2022. [15] Stefan Elfwing, Eiji Uchibe, and Kenji Doya. Sigmoid-weighted linear units for neural network function approximation in reinforcement learning. Neural networks, 107:3–11, 2018. [16] Chris Finlay, Jörn-Henrik Jacobsen, Levon Nurbekyan, and Adam Oberman. How to train your neural ode: the world of jacobian and kinetic regularization. In International conference on machine learning, pages 3154–3164. PMLR, 2020. [17] Jean Gallier. Remarks on the cayley representation of orthogonal matrices and on perturbing the diagonal of a matrix to make it invertible. arXiv preprint math/0606320, 2006. [18] Ian Gemp, Brian McWilliams, Claire Vernade, and Thore Graepel. Eigengame: Pca as a nash equilibrium. International Conference on Learning Representations, 2021. [19] Jan J Gerbrands. On the relationships between svd, klt and pca. Pattern recognition, 14(1-6): 375–381, 1981. [20] Roger G Ghanem and Pol D Spanos. Stochastic finite elements: a spectral approach. Courier Corporation, 2003. [21] Amir Gholaminejad, Kurt Keutzer, and George Biros. Anode: Unconditionally accurate memoryefficient gradients for neural odes. In Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI-19, pages 730–736. International Joint Conferences on Artificial Intelligence Organization, 7 2019. doi: 10.24963/ijcai.2019/103. URL https: //doi.org/10.24963/ijcai.2019/103. [22] Gene H Golub and Charles F Van Loan. Matrix computations. JHU press, 2013. [23] David Gottlieb and Steven A Orszag. Numerical analysis of spectral methods: theory and applications. SIAM, 1977. [24] Brian C Hall. Lie groups, lie algebras, and representations. In Quantum Theory for Mathematicians, pages 333–366. Springer, 2013. [25] Nicholas J Higham. Accuracy and stability of numerical algorithms. SIAM, 2002. [26] Kurt Hornik, Maxwell Stinchcombe, and Halbert White. Multilayer feedforward networks are universal approximators. Neural networks, 2(5):359–366, 1989. [27] Arieh Iserles. On cayley-transform methods for the discretization of lie-group equations. Foundations of Computational Mathematics, 1(2):129–160, 2001. [28] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014. [29] Andreas M Kist. CelebA Dataset cropped with Haar-Cascade face detector, 2021. URL https://doi.org/10.5281/zenodo.5561092. [30] Bernard O Koopman. Hamiltonian systems and transformation in hilbert space. Proceedings of the National Academy of Sciences, 17(5):315–318, 1931. [31] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. Journal of Machine Learning Research, 24(89):1–97, 2023. [32] Alex Krizhevsky and Geoffrey Hinton. Learning multiple layers of features from tiny images. Technical Report, 2009. [33] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278–2324, 1998. 11

[34] Bruno Lévy. Laplace-beltrami eigenfunctions towards an algorithm that" understands" geometry. In IEEE International Conference on Shape Modeling and Applications 2006 (SMI’06), pages 13–13. IEEE, 2006. [35] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations. arXiv preprint arXiv:2003.03485, 2020. [36] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. International Conference on Learning Representations, 2021. [37] Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar. Physics-informed neural operator for learning partial differential equations. ACM/JMS Journal of Data Science, 1(3):1–27, 2024. [38] Lu Lu, Pengzhan Jin, Guofeng Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 3(3):218–229, 2021. [39] Stéphane Mallat. A wavelet tour of signal processing. Elsevier, 1999. [40] Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng. Nerf: Representing scenes as neural radiance fields for view synthesis. Communications of the ACM, 65(1):99–106, 2021. [41] Vismay Modi, Nicholas Sharp, Or Perel, Shinjiro Sueda, and David IW Levin. Simplicits: Mesh-free, geometry-agnostic elastic simulation. ACM Transactions on Graphics (TOG), 43(4): 1–11, 2024. [42] Aviv Navon, Aviv Shamsian, Idan Achituve, Ethan Fetaya, Gal Chechik, and Haggai Maron. Equivariant architectures for learning in deep weight spaces. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 25790–25816. PMLR, 23–29 Jul 2023. URL https://proceedings.mlr.press/v202/navon23a.html. [43] David Pfau, Stig Petersen, Ashish Agarwal, David GT Barrett, and Kimberly L Stachenfeld. Spectral inference networks: Unifying deep and spectral learning. International Conference on Learning Representations, 2019. [44] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686–707, 2019. [45] Yousef Saad. Numerical methods for large eigenvalue problems: revised edition. SIAM, 2011. [46] Peter J Schmid. Dynamic mode decomposition of numerical and experimental data. Journal of fluid mechanics, 656:5–28, 2010. [47] Nicholas Sharp, Cristian Romero, Alec Jacobson, Etienne Vouga, Paul Kry, David IW Levin, and Justin Solomon. Data-free learning of reduced-order kinematics. In ACM SIGGRAPH 2023 Conference Proceedings, pages 1–9, 2023. [48] Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Implicit neural representations with periodic activation functions. Advances in neural information processing systems, 33:7462–7473, 2020. [49] Henry Stark and John W Woods. Probability, random processes, and estimation theory for engineers. Prentice-Hall, Inc., 1986. [50] Petre Stoica, Randolph L Moses, et al. Spectral analysis of signals, volume 452. Pearson Prentice Hall Upper Saddle River, NJ, 2005. 12

[51] Matthew Tancik, Pratul Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan Barron, and Ren Ng. Fourier features let networks learn high frequency functions in low dimensional domains. Advances in neural information processing systems, 33:7537–7547, 2020. [52] Geoffrey Ingram Taylor and Albert Edward Green. Mechanism of the production of small eddies from large ones. Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences, 158(895):499–521, 1937. [53] Lloyd N Trefethen. Spectral methods in MATLAB. SIAM, 2000. [54] Jane-Ling Wang, Jeng-Min Chiou, and Hans-Georg Müller. Functional data analysis. Annual Review of Statistics and its application, 3:257–295, 2016. [55] Yiheng Xie, Towaki Takikawa, Shunsuke Saito, Or Litany, Shiqin Yan, Numair Khan, Federico Tombari, James Tompkin, Vincent Sitzmann, and Srinath Sridhar. Neural fields in visual computing and beyond. In Computer graphics forum, volume 41, pages 641–676. Wiley Online Library, 2022. [56] Dejia Xu, Peihao Wang, Yifan Jiang, Zhiwen Fan, and Zhangyang Wang. Signal processing for implicit neural representations. Advances in Neural Information Processing Systems, 35: 13404–13418, 2022. [57] Gizem Yüce, Guillermo Ortiz-Jiménez, Beril Besbinar, and Pascal Frossard. A structured dictionary perspective on implicit neural representations. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 19228–19238, 2022. [58] Biao Zhang and Rico Sennrich. Root mean square layer normalization. Advances in neural information processing systems, 32, 2019.

Proof of Theorem 1

We first restate the theorem for completeness: Theorem 1 (Approximating the Orthogonal Group via Rank-2 Generators). Let QT denote the set of all operators Q(T ) obtained by solving the initial value problem (5) at time T > 0, where the generator K(t) is uniformly Hilbert-Schmidt, continuous, and rank-2. Then, QT is dense in SO(HS) with respect to the Hilbert-Schmidt norm topology. In other words, every rotation in this infinite orthogonal group is arbitrarily well-approximated by some time-evolving rank-2 generator. Here, we repeatedly use the fact that any rank-2 operator K(t) can be written in the following form: K(t) = α(t) ⊗ β(t) − β(t) ⊗ α(t),

(16)

where α, β : [0, T ] → H are H-valued functions of time. Our goal is to ultimately show that for any operator Q ∈ SO(HS), there exist α, β such that the solution to the following initial value problem closely approximates Q at t = T : Q̇(t) = [α(t) ⊗ β(t) − β(t) ⊗ α(t)] Q(t),

Q(0) = I.

(17)

To do so, we first restrict our space to the following space of operators: Definition 2. Let SO(fin) ⊂ SO(HS) denote the group of “finite-like” orthogonal Hilbert Schmidt operators that act as identity outside of some finite-dimensional linear subspace of V of H. That is, if Q ∈ SO(fin) then there exists some finite n and a finite set S = {e1 , . . . , en } ⊂ H such that Qv = v

for all v ∈ V ⊥ = span{e1 , . . . , en }⊥ .

We now constructively prove that the approximation theorem holds for SO(fin): Lemma 3. For any Q ∈ SO(fin), there exists continuous maps α, β : [0, T ] × Ω → R such that the solution Q(T ) of (17) satisfies Q(T ) = Q. Proof. Let Q ∈ SO(fin). Then, by definition, there exists a finite-dimensional subspace V alongside its orthonormal basis {e1 , . . . , en } ⊂ H such that Q acts trivially on V ⊥ . We can therefore restrict the representation of Q to V and get a finite representation QV ∈ SO(n). By classical results in numerical linear algebra and Lie group theory [24], any element of SO(n) can be written as a finite product of Givens’ rotations acting on coordinate planes. In other words, for any QV there exist M angles {θm }M m=1 and pairwise orthogonal directions M n−1 {am }M such that: m=1 , {bm }m=1 ⊂ S QV =

M Y

RaGivens (θm ) = m ,bm

m=1

⊤ exp θm (am b⊤ m − bm am ) .

(18)

m=1

√ √ Pn Pn We can now define αm = θm · i=1 ami ei ∈ H and βm = θm · i=1 bmi ei ∈ H which allows us to lift the representation back from the coefficient space to H and write: Q=

M Y

exp αm ⊗ βm − βm ⊗ αm .

(19)

m=1

Now, the idea is to construct a continuous field in H, that once integrated, results in the operator above. We achieve this by stitching together M time-evolving functions below: αm (t) = rm (t) αm ,

βm (t) = rm (t) βm .

(20)

We define rm : R → R as a continuous control function in the form of a “bump” that has the R1 2 following properties: rm (0) = rm (1) = 0 and 0 rm (τ )dτ = 1. One such control function can be the square root of the probability density of the Beta distribution. These properties on rm (t) ensure αm (1) = αm+1 (0) and βm (1) = βm+1 (0), allowing us to stitch generators continuously and form a bigger generator. Furthermore, the construction above allows the self-adjoint operator to always self-commute over all t ∈ [0, 1]: 2 αm (t) ⊗ βm (t) − βm (t) ⊗ αm (t) = rm (t) [αm ⊗ βm − βm ⊗ αm ]

(21)

Next, w.l.o.g. assume T = M and define the following M initial value problems: Qm−1 (1) if m > 1 Q̇m = [αm (t) ⊗ βm (t) − βm (t) ⊗ αm (t)] Qm (t), Qm (0) = . I otherwise

(22)

Using the fact that the generator above is skew-adjoint and self-commutes, we can write the solution of the ODE as an exponential map and simplify:   Z 1     2 Qm (1) = exp [αm ⊗ βm − βm ⊗ αm ] · rm (τ )dτ · Qm (0)   | 0 {z }

(23)

(=1) by construction

⇒ QM (1) =

M Y

exp (αm ⊗ βm − βm ⊗ αm )

(24)

m=1

Hence, chaining them yields a time-evolving Q(t) for the form (17) where Q(T ) = QM (1) = Q. Lemma 4. SO(fin) is dense in SO(HS) with respect to the Hilbert-Schmidt operator norm. Proof. Let Q ∈ SO(HS), then from the definition, it can be written as the exponential map of a Hilbert-Schmidt skew-adjoint operator, i.e., Q = exp(K) for some K = −K ∗ where ∥K∥HS < ∞. Since K is Hilbert-Schmidt, using the spectral theorem for operators, there exists a sequence of real values λn > 0 and an orthonormal set of vectors in H, {un , vn }∞ n=1 , such that: K=

∞ X

λj (vj ⊗ uj − uj ⊗ vj ).

(25)

j=1

Next, we can define the sequence of finite rank operators {Kn }∞ n=1 as the spectral truncation: Kn =

n X

λj (vj ⊗ uj − uj ⊗ vj ).

(26)

j=1

Importantly, limn→∞ Kn = K in ∥ · ∥HS and K commutes with all Kn : KKn = Kn K. Since Kn is finite rank, its exponential acts as identity on range(Kn )⊥ ; this means that Qn := exp(Kn ) ∈ SO(fin). Furthermore, we have that: lim ∥Qn − Q∥HS

(27)

n→∞

= lim ∥ exp(Kn ) − exp(K)∥HS

(28)

= lim ∥ exp(K) [exp(Kn − K) − I] ∥HS Exponentiation of Commutative Operators n→∞ ≤ ∥ exp(K)∥HS · ∥ exp lim Kn − K − I∥HS Continuity of norms and exponentials

(29)

= ∥ exp(K)∥HS · ∥ exp(0) − I∥HS = 0.

(31)

n→∞

(30)

Thus, we have found for any Q ∈ SO(HS), a sequence of operators Qn ∈ SO(fin) such that n→∞ Qn =⇒ Q in ∥ · ∥HS ; this is precisely the definition of SO(fin) being dense in SO(HS). Collectively, the two lemmas above prove Theorem 1. Finally, a universal approximation result with MLPs is straightforward: with a sufficiently overparameterized θ, any two MLPs can approximate both functions α and β on the compact domain [0, T ] × Ω [26]. In turn, by approximating the continuous α : [0, T ] × Ω → R and β : [0, T ] × Ω → R, these MLPs implicitly approximate an arbitrary Q ∈ SO(HS).

Diagonalizing Linear Operators

B.1

Proof of Theorem 2

∞ Theorem 2. Let A be a trace-class positive semidefinite linear i=1 be a P∞ operator, and let {p(i)} ∞ strictly decreasing sequence of positive weights such that i=1 p(i) = 1, and {φi }i=1 be any orthonormal basis on H. Consider Qopt as the solution to the following variational problem: X max p(i) · ⟨Qφi , AQφi ⟩ (7) Q∈SO(HS) i=1

Then, ϕ⋆i := Qopt φi are the eigenfunctions of A. That is, Aϕ⋆i = λi ϕ⋆i for all i ≥ 1, where λ1 ≥ λ2 ≥ λ3 ≥ · · · are the corresponding eigenvalues. Proof. Let Υ and Λ be diagonal operators with entries p(i) and λi , respectively, and let Φ ∈ SO(∞) be the orthogonal operator whose i-th column corresponds to Qφi . We can compactly rewrite the objective as tr(ΥΦ∗ AΦ). Let σi (·) an operator functional that inputs an operator and returns its ith largest eigenvalue. Since both Υ and A are trace-class, by applying von Neumann’s trace inequality and the Poincaré separation theorem, respectively, we upper-bound the objective as follows: ∞ X tr(ΥΦ∗ AΦ) ≤ p(i)σi (Φ∗ AΦ) Von-Neumann inequality for trace-class operators (32) i=1

≤ =

∞ X i=1 ∞ X

Φ is unitary, so: σi (Φ∗ AΦ) ≤ σi (A)

p(i)σi (A) p(i)λi .

(33) (34)

i=1

The upper bound is exactly achieved when Φ ∈ SO(HS) diagonalizes A. B.2

Diagonalizing the Covariance Operator

An important detail omitted from the main text is that, in practice, many datasets X ∼ P used for principal component analysis are not necessarily zero-mean. To properly diagonalize the true covariance operator C using Qθ , we must shift the data by its empirical mean. Rather than computing the empirical mean beforehand, we accomplish the same effect by concurrently fitting a spatial mean function—parameterized as a neural field νψ : Ω → R using random Fourier features [51]. This leads to a joint objective where we maximize the captured energy of the centered data while minimizing the L2 reconstruction error of the mean: max Ei∼p,X∼P ⟨X − νψ , Qθ φi ⟩2 , (35) θ

min EX∼P ∥νψ − X∥2 . ψ

(36)

Because the parameters ψ and θ are decoupled across these two objectives, we can optimize them in parallel using independent optimizers for each task. For clarity, the algorithm is summarized as a pseudo-code in algorithm 1. B.3

Diagonalizing the Neural Tangent Kernel

Let MLP(·; ξ) : Ω → R a non-linear neural function parameterized by ξ. We restrict our attention to networks with a scalar output; this encompasses the experimental setting in section 7 like our Two Moons experiment where the output of the network represents the logit of a binary classifier. The Neural Tangent Kernel (NTK) is defined by the inner product of the network’s parameter gradients. We can express this as an integral operator κξ acting on the function space H with the kernel κξ [x, y] = ∇ξ MLP(x; ξ)⊤ ∇ξ MLP(y; ξ). Specifically, the action of the NTK operator on any function f ∈ H is defined as: Z (κξ f )(x) = ∇ξ MLP(x; ξ)⊤ ∇ξ MLP(y; ξ)f (y) µ(dy). (37) Ω

Algorithm 1: Continuous Principal Component Analysis Input: Dataset X , reference basis {φi }i=1 , distribution p(i), domain Ω with measure µ. Output: Optimized basis Qθ and mean function νψ . Initialize orthogonal network parameters θ and mean function parameters ψ; while not converged do b ∼ µ; Sample a function X ∼ X , a basis index i ∼ p, and quadrature points Ω P 1 2 Lmean ← |Ω| b ∥νψ (ω) − X(ω)∥ ; b ω∈Ω P 1 2 JPCA ← |Ω| b ⟨X(ω) − νψ (ω), (Qθ φi )(ω)⟩ ; b ω∈Ω ψ ← update(ψ, ∇ξ Lmean ); θ ← update(θ, ∇θ JPCA );

To diagonalize the NTK using our framework, we substitute κξ for A in the variational objective of Theorem 2. For a basis element ϕi = Qφi , we can simplify the quadratic form ⟨ϕi , κξ ϕi ⟩: ⟨ϕi , κξ ϕi ⟩ Z Z ⊤ = ϕi (x) ∇ξ MLP(x; ξ) ∇ξ MLP(y; ξ)ϕi (y) µ(dy) µ(dx) Ω

(38) (39)

Ω

⊤ Z

[∇ξ MLP(x; ξ) · ϕi (x)] µ(dx)

[ϕi (x) · ∇ξ MLP(x; ξ)] µ(dx)

= Ω

(Fubini)

(40)

Ω

= ∥ Ex∼µ [ϕi (x)∇ξ MLP(x; ξ)] ∥22 .

(41)

Thus, the quadratic form is in fact equal to the squared norm of the expected parameter gradient. Substituting this simplified quadratic form back into the general objective, we arrive at our specialized objective function for diagonalizing the NTK: 2

max Ei∼p ∥Ex∼µ [(Qθ φi )(x)∇ξ MLP(x; ξ)]∥2 . θ

(42)

In practice, the optimization proceeds by sampling a target index i ∼ p and a batch of quadrature b ∼ µ. At these points, we evaluate both our trial basis function (Qθ φi )(x) and the per-sample points Ω parameter gradients ∇ξ MLP(x; ξ). We then compute the empirical expectation over the batch by taking the weighted average of these evaluated gradients, yielding a single vector in the network’s parameter space. Finally, we compute the squared norm of this expectation and backpropagate to update the basis parameters θ, keeping the MLP parameters ξ frozen. Upon convergence, the learned functions ϕi = Qθ φi faithfully represent the dominant eigenmodes of the training dynamics. For clarity, the algorithm is summarized as a pseudo-code in algorithm 2. Algorithm 2: Diagonalizing the Neural Tangent Kernel Input: Domain Ω with measure µ, reference basis {φi }i=1 , distribution p(i), frozen MLP(·; ξ). Output: Optimized NTK eigenfunctions ϕi = Qθ φi . Initialize orthogonal network parameters θ; while not converged do b ∼ µ and a basis index i ∼ p; Sample a batch of quadrature points Ω // Step 1: Evaluate trial eigenfunctions and parameter gradients b do for each ω ∈ Ω g(ω) ← ∇ξ MLP(ω; ξ); q(ω) ← (Qθ φi )(ω); // Step 2: Compute empirical expected gradient and NTK objective P 1 v ← |Ω| b q(ω)g(ω); b ω∈Ω JNTK ← ∥v∥22 ; // Step 3: Update basis generator θ ← update(θ, ∇θ JNTK );

Cayley Integration Details and Multichannel Details

This appendix expands on section 6. We provide a pseudo-code for the Cayley integration in b ⊂ Ω. algorithm 3; this method is internally implemented whenever Qθ φi (ω) is called for ω ∈ Ω Notably, we set T = 1 for all experiments. We also extend the Cayley integration to multichannel function spaces and provide more details on the forward and backward Euler ODE solvers used in our calculations in Figure 2. Algorithm 3: Continuous Change-of-Basis via Cayley Integration Input: Initial basis function φi , parameters θ, rank r. Output: Estimated transformed basis evaluation ϕb ≈ Qθ φi . Sample random time discretization 0 = t0 < t1 < · · · < tL = T ; b = {ωj }D ∼ µ; Sample random dense spatial samples Ω j=1 b Initialize ϕb ← φi (Ω); for ℓ = 1, . . . , L do ∆tℓ ← tℓ − tℓ−1 and tmid ← tℓ−1 + ∆tℓ /2; √ b ← ∆tℓ Uθ (tmid , Ω) b ∈ Rr×D and Sb ← Mθ (tmid ) − Mθ (tmid )⊤ ∈ Rr×r ; U // Step 1: Estimate r × r Gram matrix and Woodbury inverse M −1 b i (ω)U b j (ω) ∈ Rr×r and M ← Ir×r − 1 G Sb bU b ⊤ ]ij ≈ 1 P b U G ← [U ; ω∈Ω D 2 1 b⊤ b b b b b // Step 2: Apply Forward Action ϕ 7→ ϕ + 2 U S U ϕ r b b ϕb ≈ 1 P b U b c←U ω∈Ω (ω)ϕ(ω) ∈ R ; D 1 ⊤ b Sb c; yb ← ϕb + 2 U // Step 3: Apply Backward Action via Woodbury b yb ≈ 1 P b U b c̃ ← U y (ω) ∈ Rr ; ω∈Ω (ω)b D b c̃ ∈ Rr ; w ← SM b b ⊤ w; ϕ ← yb + 12 U C.1

Multichannel Hilbert Spaces

The multichannel-valued Hilbert space HC = L2 (Ω, RC ), is defined through the following inner product that integrates over the spatial domain and sums across channels: C Z X ⟨f, g⟩ = fc (x)gc (x) dµ(x). (43) c=1

Ω

Consequently, HC is a single unified Hilbert space rather than C independent ones. Naturally, when C = 1, HC simplifies to the standard real-valued Hilbert space denoted as H in the main text. For experiments involving multichannel data (e.g. RGB images in CelebA), we only need a minor adjustment to Uθ . Instead of the single-channel formulation Uθ (t, ω) ∈ Rr , we use Uθ (t, ω) ∈ RrC bci if it is evaluated and index by Uθ (t, ω)ic to denote the cth channel of the ith row of Uθ , or in short U b Inner product approximations are updated to account for the at some time and quadrature points Ω. channel dimension, meaning the Gram matrix and matrix-vector products become: h i 1 X b bU b⊤ b (ω)jc , U (ω)ic · U (44) U ≈ b ij |Ω| b ω∈Ω 1≤c≤C

i X b fb ≈ 1 b (ω)ic fc (ω). U U b i |Ω| b

(45)

ω∈Ω 1≤c≤C

If, instead, we applied generators independently per channel, we would have forced a block-diagonal structure that isolated the color channels. This would have prevented us from transferring energy 18

Solver exploded Forward Euler Backward Euler Cayley Transform

Backward Euler energy decreases

Cayley Transform Forward energy Euler is preserved

energy increases

Constant energy sphere

(a) Forward integration of constant Fourier basis

(b) Gram matrix during evolution

Figure 7: Full Cayley Ablation. Comparison of the Cayley method against backward and forward Euler. (a) Example evolution of a basis function ϕθi (t) over time. The Cayley solver perfectly maintains signal energy, while forward Euler rapidly explodes and backward Euler dissipates. (b) Beyond norm preservation, the inner products of various basis functions, visualized as a gram matrix, also remain preserved through various snapshots of time. (c) A conceptual diagram of the Lie group integrator, which strictly remains on the orthogonal manifold, compared to methods that either overshoot or decay. Overall, Cayley integration is critical for stable optimization.

between channels, instead, by operating on the unified multichannel Hilbert space HC , our formulation allows mixing between channels. Empirically, this is what allows us to map a blue Fourier basis element to a multicolored function, as shown in Figure 2. Moreover, since HC itself is a proper Hilbert space, the approximation Theorem 1 also holds in the multichannel case; our entire proof transfers by replacing the definition of the inner product. C.2

Cayley vs. Off-the-Shelf ODE Solvers

Figure 7 compares our Cayley integrator with two classical alternatives, holding the number of steps b ℓ constant: L and the skew-symmetric operator K b ℓ )ϕbℓ−1 . For • Forward (explicit) Euler: The update formula can be written as ϕbℓ = (Ip+ ∆t K b a skew-symmetric Kℓ , this is unstable as each step amplifies the norm by 1 + (∆t σK )2 > 1, b ℓ . This leads to the explosive growth shown in where σK represents the singular values of K Figure 2 and Figure 7. b ℓ )−1 ϕbℓ−1 . This • Backward (implicit) Euler: The update rule can be written asp ϕbℓ = (I − ∆t K approach is stable but dissipative, with the norm decaying by 1/ 1 + (∆t σK )2 < 1 at each step. This leads to the dissipative performance shown in Figure 2 and Figure 7. For a fair comparison, we implement the inverse using the same Woodbury identity as our Cayley solver. • Cayley (ours): Norm-preserving, ensuring ∥ϕbℓ ∥ = 1 at every integration step.

Experimental Details

This appendix collects the configurations used for every experiment in section 7. All runs share the same Cayley integrator (Appendix C) and differ only in their variational objective (PCA, NTK diagonalisation, or Koopman fitting). Table 1 also summarized some key hyperparameters that are different across the different experiments. D.1

Indexing the Infinite Bases and their Prior p(i)

All variational problems in section 5 involve a real Fourier basis φi defined on Ω = [0, 1]d and an index prior p(i). While Fourier bases are enumerated by integers, for higher-dimensional domains, they are often enumerated with vector multi-indices. When d > 1, the index i is represented as a vector of size d, i = [i1 , i2 , . . . , id ]. For multichannel cases, the index is additionally extended to size d + 1, i = [i1 , i2 , . . . , id , c], for multi-channel tasks. For each spatial dimension j, the absolute value |ij | denotes the Fourier frequency. The sign of ij determines the trigonometric function: ij > 0 yields a cosine, ij < 0 yields a sine, and ij = 0 yields the constant function. For example, in a 2D single-channel setting, the index i = [2, −1] corresponds to the scalar basis function φi (x1 , x2 ) ∝ cos(2 × 2πx1 ) sin(1 × 2πx2 ). In a multi-channel setting (e.g., a 3-channel RGB image), the extended index i = [2, −1, 3] applies this exact same spatial feature strictly to the third channel. This yields the vector-valued basis function φi (x1 , x2 ) ∝ [0, 0, cos(2 × 2πx1 ) sin(1 × 2πx2 )]. We quantify the overall spatial frequency of a Fourier basis element using the squared ℓ2 -norm of its Pd spatial multi-index i = [i1 , . . . , id ], given by ∥i∥22 = j=1 i2j . To prioritize lower-frequency features, we define the index prior p as a strictly decreasing power law over this norm: p(i1 , . . . , id , c) =

−α 1 1 + ∥i∥22 , C Zα

Zα :=

1 + ∥i∥22

−α

(46)

k∈Zd

where C is the total channel count and α > d/2 ensures the partition function Zα is finite. Because Theorem 2 requires a strict total ordering, any ties between distinct basis elements with identical frequency norms (e.g., i = (2, 1) and i = (1, 2)) are resolved by applying a deterministic tiebreaker. Objective Variance Reduction. Naïve Monte-Carlo estimation of variational objectives involving the expectation Ei∼p [·] can suffer from high variance. To mitigate this, we use a stratified-with-tail estimator that decomposes the expectation into two parts: an exact stratum S = {i : p(i) ≥ τ } comprising all high-probability indices, and a Monte-Carlo tail consisting of Nτ independent and identically distributed (i.i.d.) draws from p, explicitly rejecting any indices already present in the exact stratum. Concretely, for an integrand f , the estimator is N

b ]= E[f

p(i) f (i) +

i∈S

exact stratum

}

τ 1 X f (i(k) ) 1[i(k) ∈ / S], Nτ k=1 | {z }

i.i.d.

i(k) ∼ p,

(47)

Monte-Carlo tail

where 1(·) is the indicator function of an event. By construction, this estimator is unbiased: the exact stratum is evaluated with zero variance, while the variance of the sampled tail scales as O(1/Nτ ). The probability threshold τ and the tail sample size Nτ for each experiment are detailed in Table 1. D.2

Architecture

All of our experiments parameterize the finite-rank basis generator using two neural networks: Mθ (t), which controls the internal mixing of the basis vectors, and Uθ (t, ω), which projects the continuous spatial domain into the rank-r subspace. Both networks are conditioned on time via a shared sinusoidal positional embedding γ(t) ∈ R512 . This embedding consists of F = 256 frequencies log-spaced in [20 , 23 ]: γ(t) = [sin(2πω1 t), cos(2πω1 t), . . . , sin(2πωF t), cos(2πωF t)] , 20

ωf = 23(f −1)/(F −1) .

(48)

Table 1: Hyperparameters for every experiment shown in section 7. Rank r, time-discretisation L, and quadrature size D are the parameters that control accuracy of the function-space ODE; τ is the tail cutoff and Nτ the tail size for the stratified energy estimator. Experiment 1-D PCA (Figure 3) INR-Zoo MNIST (Figure 4) CelebA (Figure 4) NTK two-moons (Figure 5) Fluid Vortices (Figure 6)

r 30 20 50 10 10

L 20 20 20 20 50

Nτ

64 642 642 322 322

10−3

16 64 32 32 32

10−4 10−5 7×10−4 10−4

Internal Mixing Matrix Mθ (t). The network Mθ : R512 → Rr×r maps the time embedding γ(t) to the unconstrained mixing matrix. We use a small two-hidden-layer MLP with hidden widths [64, 64], SiLU activations [15], and RMSNorm [58] applied between the linear layers and activations. Because only the skew-symmetric component Mθ (t) − Mθ (t)⊤ enters the dynamics ODE, the raw output of this network does not require any structural constraints. Spatial Projection Field Uθ (t, ω). The network mapping the spatio-temporal inputs (γ(t), ω) 7→ U ∈ Rr dictates how the spatial domain evolves. To efficiently capture high-frequency details while ensuring robust frequency mixing and representing the 1/k 2 power spectrum typical of natural signals, we construct the final projection field Uθ (t, ω) as the sum of a base NeRF-style multi-layer perceptron [40] UθMLP (t, ω) and a time-varying basis-element-selector residual: X 1 φk (ω). (49) Uθ (t, ω) = UθMLP (t, ω) + wk,θ (t) 1 + ∥k∥22 |k|∞ ≤K−1

Across all experiments, the base network UθMLP (t, ω) concatenates [γ(t), ω, η(ω)], where η(ω) represents multi-scale random Fourier features. We use L = 8 levels log-spaced over frequencies [1, 64], with nℓ = max(1, 64/2ℓ ) projections per level (yielding ≈ 254 features). This input is fed into a 4-layer MLP of width 256 with SiLU activations and RMSNorm, leading to a linear map to Rr . For the residual component, the time-varying coefficients wk,θ : R512 → Rr are parameterized by a small 2-layer MLP operating on the time embedding γ(t) (hidden widths [64, 64], SiLU, RMSNorm). The 1/(1 + ∥k∥22 ) factor acts as a structural low-pass prior, while the bandwidth hyperparameter K is set on a per-experiment basis. D.3

Compute and Optimization

All experiments run on a single NVIDIA RTX 5090 (32 GB) using the Adam optimizer [28] with default hyperparameters and learning rate 10−3 , gradient clipping at ∥∇∥2 ≤ 1, and gradient checkpointing through the Cayley step ladder to keep VRAM constant in L. Each Cayley step costs O(r3 + r2 D); for the typical 1024 ≤ D ≤ 4096, 10 ≤ r ≤ 50 this is well below the cost of a single b MLP forward pass through Uθ , so wall-clock is dominated by the L network evaluations of Uθ (tℓ , Ω) per training step. A typical 2-D run (CelebA RGB) trains in ∼ 1.5 hours. D.4

Additional Experiments

Additional experiments and visualizations are provided in Figure 8–Figure 11. Figure 8 evaluates the learned basis on the CIFAR10 SIREN Zoo [32], showing that the learned bases resemble smooth principal components of the dataset and yield higher PSNRs and greater captured energy than the Fourier basis. Figure 9 presents analogous results for CelebA, where the learned bases take the form of smooth eigenfaces and achieve better reconstruction quality and higher cumulative captured energy with fewer basis elements. Figure 10 visualizes a truncated set of initial and learned bases across applications—including MNIST, Taylor-Green vortices, and CelebA—demonstrating that while the initialization corresponds to a random path on the Lie manifold, training recovers the eigenmodes of the covariance operator. Finally, Figure 11 illustrates the internal network dynamics for the vortex experiment, visualizing the learned neural fields Uθ0...4 and spatial rotations exp(Mθ − Mθ⊤ ) that define the generator Kθ (t) and enable simulation of the Koopman operator via rank-r planar rotations of the Fourier basis. 21

Fourier Learned basis (Ours)

Cumulative captured energy

Learned basis (Ours)

Fourier

(a) CIFAR10 INR (SIREN zoo) dataset

Fourier Learned (Ours)

number of basis elements (c) CIFAR10 INR cumulative captured energy

(b) CIFAR10 INR reconstruction

Figure 8: PCA for CIFAR10 [32] SIREN Zoo. (a) Visualizes the reference and learned basis after training; the learned bases resemble smooth principal components of the CIFAR10 dataset. (b) Reconstruction with the learned basis results in higher PSNRs. (c) Quantitatively, the captured energy is higher for our learned basis compared to the Fourier basis.

Fourier Learned basis (Ours) Fourier Learned basis (Ours)

Cumulative captured energy

Learned basis (Ours)

Fourier

(a) CelebA dataset

Fourier Learned (Ours)

number of basis elements (c) CelebA cumulative captured energy

(b) CelebA reconstruction

Figure 9: PCA for CelebA. (a) Visualizes the reference and learned basis; the learned bases resemble smooth “eigenfaces” of CelebA. (b) Reconstruction with the learned basis obtains higher PSNRs and captures the structure of the face better with fewer elements. (c) Quantitatively, the cumulative captured energy is higher for our learned basis compared to Fourier basis.

Fourier basis

Learned basis (Ours)

(a) MNIST INR dataset covariance operator learning basis mosaic

Initialized basis for learning

Learned basis (Ours)

Fourier basis

(b) Taylor-Green Koopman operator learning basis mosaic

Learned basis (Ours)

Figure 10: Initial and Learned 2D and Multichannel Bases. Visualization a truncated set of initial and learned bases across applications using our index notation (see subsection D.1); the center of the mosaic corresponds to the direct current (DC) Fourier functions. (a) 2D Fourier and a learned bases fit to the principal components of the MNIST SIREN Zoo; the bases are visualized as a 2D mosaic indexed as φi1 ,i2 for Fourier and Qθ φi1 ,i2 for the learned basis where i1 and i2 correspond to the horizontal and vertical frequencies. (b) Similar results for the Vortices where each Qθ φi1 ,i2 corresponds to a dynamic mode of the Taylor-Green vortex field. (c) The 2D RGB-channel Fourier and learned basis for CelebA. Here the bases are visualized as a 2D mosaic with an extra channel dimension and indexed as φi1 ,i2 ,c . The basis Qinit θ φi1 ,i2 ,c at initialization is obtained from a random path on the Lie manifold, but once trained, it learns to map to the eigenmodes of the covariance operator.

(rows: ranks, columns: time)

Row

(a) Learned generators

Column

Column (b) Rank rotation matrices

Column

Figure 11: Internal Network Dynamics. This figure illustrates the learned neural network outputs that transform the Fourier basis into the dynamic modes of the Taylor-Green vortices. Specifically, we visualize the first five neural fields, Uθ0...4 , across various time steps (a), alongside the spatial rotations exp(Mθ − Mθ⊤ ) (b) that jointly characterize the generator Kθ (t). Passing the Fourier basis through these rank-r planar rotations allows us to simulate the Koopman operator.

Record · ID 204802 · SHA-256 fe47aaa41643f683
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.