Score the Algebra, Not the Span: Dimension Reduction for Transfer Operator Models of Dynamical Systems

Score the Algebra, Not the Span: Dimension Reduction for Transfer Operator Models of Dynamical Systems

arXiv:2608.18918v1 [cs.LG] 19 Aug 2026

Mark Kozdoba1 Technion, IIT

Shie Mannor Technion, IIT and NVIDIA

Abstract Dimension reduction for dynamical systems is standard practice, and the standard route is spectral: model the transfer (Koopman) operator by its leading modes. We show that on systems assembled from several weakly interacting components — a structure common in physical and biological settings — this may either require an exponential number of modes, or drop an entire component: the component is absent from the model rather than modeled coarsely, and no function of it can be predicted at any accuracy. We call this linear masking. The cause is that a rank-based model pays one coordinate per mode. We propose to score instead the σ-algebra the coordinates generate, so that products and powers come free and a component’s cost is governed only by its generators rather than by all its interactions. The criterion is a χ2 -divergence between the embedded present and future, and it carries a budget guarantee: twice the intrinsic dimension of the dynamics is enough coordinates for an embedding whose algebra carries the operator’s entire spectrum, with its full infinite rank. In variational form the criterion admits off-the-shelf estimators, and restricting its critic to the bilinear class returns the VAMP score on the span, so rank-based methods are one end of the same family. We demonstrate the proposed objective on a composite of published benchmark systems. We exhibit examples where the rank-based methods completely miss the masked components at all ranks k < 100, while ten algebra coordinates recover all of them. In addition, the resulting algebra representation supports predicting the masked components from few labels, while direct regression from the high-dimensional observation or from the VAMP features fail.

1

Introduction

Learning a model of a dynamical system from observed trajectories is a basic machine learning task. Instead of modeling the dynamics on the state space X directly, one standard approach is to work with the linear operator T that the dynamics induces on the space of functions f over X , carrying a function of the future state to its conditional expectation given the present. In deterministic systems this is known as the Koopman operator [Brunton et al., 2022], and it is simply the transition operator in the Markov setting [Meyn and Tweedie, 2009]. For many physical systems the eigenfunctions of T are known in closed form, and serve as the basis of the analysis of the dynamics. However, when the system is given only through sampled trajectories, the eigenfunctions have to be learned, and the standard target to be approximated is a finite-dimensional invariant subspace. For instance, the extended dynamic mode decomposition (EDMD) fits one inside a chosen dictionary [Williams et al., 2015a], while in the stochastic and non-reversible setting, one of the most commonly used objectives for learning the leading k-dimensional singular subspace is VAMP [Wu and Noé, 2020]. Its neural form, VAMPnets [Mardt et al., 2018], learns k functions that realize the best rank-k approximation to T . 1 markk@technion.ac.il

1

Multicomponent systems and linear masking. Many systems of interest are multicomponent by nature, and biological ones especially so. For instance, a biomolecular complex is well known to decompose into domains whose internal kinetics are fast and whose mutual coupling is weak. This near-product structure is then put to use: decomposing and modeling such systems domain by domain rather than jointly is an established strategy [Hempel et al., 2021, Mardt et al., 2022]. However, finding such a decomposition is itself hard. It presupposes both the number of components and the number of states allotted to each, and it applies only where those components are uncoupled or weakly coupled, a restriction that is stated explicitly in the literature [Mardt et al., 2022]. Here we ask instead whether such systems can be modeled directly, without the need to decompose the system. As mentioned above, a common approach to direct modeling of dynamical systems is by the rank-k approximation of the transfer operator T . However, the product structure above may severely limit the approximability of a system by a low-rank operator. The reason is a counting one: such a model spends one coordinate on each singular function it retains. When the components are independent, the singular functions of the joint operator are products of the singular functions of the components, and similarly the singular values are products of the components’ singular values. This implies in turn that the number of modes above a given level may grow exponentially in the number of components. We note that the most problematic manifestation of this phenomenon happens when there are spectral gaps between the components themselves. For instance, if there are a few slow components (large singular values), and a faster component (smaller leading singular value), then many products of the slow components’ values will be higher than the leading value of the fast component. In Section 5 we provide an example of a system with eight components, in which a rank k = 100 representation contains no mode involving six of them, and therefore carries no information about those six at all — they are not approximated coarsely, they are absent, while each of them is individually predictable. This is striking, since the whole system is intrinsically ten-dimensional, and one might expect a 100-dimensional representation to carry all the information about it. We call this phenomenon linear masking. Approximating by an algebra. The underlying reason for linear masking is that rank based models must assign a separate coordinate to every nonlinear interaction between existing singular functions (such as their products), even though such an interaction is already computable from the coordinates the model holds. We therefore ask for an approximation that is credited for everything computable from its coordinates, rather than one that must buy each combination separately. We formalize this intuition by looking for an embedding ϕ = (ϕ1 , . . . , ϕm ) that generates the largest σ-algebra, as measured by the amount of the spectrum of T it captures. That is, consider the space Vϕ of functions of x computable from the representation ϕ(x) alone. Then, roughly speaking, we are interested in the ϕ for which the restriction of T to Vϕ has the largest Hilbert–Schmidt norm. Although ϕ is finite dimensional, the space Vϕ is typically not, and it contains all the interactions between the coordinates ϕi . For comparison, replacing Vϕ by the span Uϕ = span {ϕ1 , . . . , ϕm } recovers the rank-m objective of VAMPnets. In Theorem 4.1 we show that a finite budget always suffices: the number of coordinates whose algebra carries the whole of T is bounded by twice the intrinsic dimension of the dynamics, a quantity that agrees with the number of degrees of freedom of the system in the cases of interest. The cost of a faithful representation is thus set by the degrees of freedom, rather than by the number of significant modes, which is what avoids the masking pathology described above. It is also worth noting that a wide literature reads a dimension off a knee in the spectrum of T [Noé and Nüske, 2013, Dsilva et al., 2018, von Lindheim, 2018, Talmon and Coifman, 2015]. Such spectrum-based notions of dimension can be arbitrarily large on a system with few degrees of freedom, as the example above shows. Practical modeling. While optimizing ϕ to score the algebra may a priori appear complex, the restricted Hilbert–Schmidt norm above is in fact exactly the χ2 dependence between the embedded present and future. Both it and the map ϕ attaining it can therefore be computed with off-the-shelf χ2 estimators [Sugiyama et al., 2012, Kanamori et al., 2009, Nguyen et al., 2010].

2

Experiments. We consider a ten-dimensional system with eight components assembled from published benchmark systems, some slow and some fast, and a parameter κ controlling the amount of dependence between them. Its trajectories are observed either through the native ten-dimensional state, or through a high-dimensional nonlinear warp, as is more common in practice. We demonstrate empirically the phenomena that the theory predicts. VAMPnets trained to convergence recover nothing of the six fast components at any rank below the product count of about one hundred, although each is individually predictable from the observations. This holds at every setting of κ and under both observation models. Maximizing the algebra criterion over ten coordinates, one per degree of freedom, recovers every component in all of these settings. The value of a learned dimension reduction often lies in the label budget: trained on abundant unlabeled transition pairs, a representation can support predicting the masked components from few labels. We show that the algebra objective produces representations where this is possible, while direct regression from the observations and rank-based representations both fail.

2

Related Work

The standard low-rank program scores the linear span of its coordinates, and so pays one coordinate per significant mode. The best rank-k span is the leading singular subspace of the transfer operator, which for a composition operator is the Koopman operator of the dynamics. Taking those leading functions as coordinates is the program of Koopman operator theory [Brunton et al., 2022], realized empirically by dynamic mode decomposition and its extended and kernel variants [Williams et al., 2015a,b]. In the stochastic, non-reversible setting the variational score for that subspace is VAMP [Wu and Noé, 2020], and its neural form, scoring the span of a learned k-dimensional feature map, is VAMPnets [Mardt et al., 2018]. The current neural mode learners target the same top-k subspace [Jeong et al., 2025, Kostic et al., 2024b, Deng et al., 2022, Pfau et al., 2019], as does the kernel and RKHS operator-regression line [Klus et al., 2020b, Kostic et al., 2022], where the rank-k truncation is replaced by a restriction on the RKHS norm. On a multicomponent system all of them over-count the intrinsic dimension by an unbounded factor, since the modes they enumerate are products and harmonics of a few generators. VAMPnet is our primary benchmark (Section 5). Closest in machinery, Turri et al. [2026] train a feature map through the same least-squares density-ratio functional we estimate with, with a bilinear critic whose optimal value they identify with the VAMP-2 score — the span once more, now reached variationally. Our one move is to leave that critic unrestricted, which scores the algebra the coordinates generate and turns the mode count into an intrinsic dimension (Section 4.2). Bittracher et al. [2018] likewise seek a mapping that preserves the dominant spectral subspace, and reach it through the embedding machinery our Theorem 4.1 also uses, over the same χ2 geometry of predictive laws. Their guarantee is a dominant-mode one by construction: the recovered eigenfunctions carry an error bounded by ε/|λi |, which diverges as |λi | → 0, so it is silent on the fast modes a rank budget discards, and they assume reversibility, which we do not. Computational mechanics reaches for the same object conceptually: causal states are the classes of a predictive-equivalence σ-algebra [Shalizi and Crutchfield, 2001], constructed in practice through kernel mean embeddings of the predictive laws [Brodu and Crutchfield, 2020], and so anchored to the kernel’s view of the observed state. Additional literature notes, on component identification and on autoencoder methods, are in Section C.

3

Preliminaries

The coupling. We work with a single state space X carrying a probability measure µ together with a Markov kernel p(· | x), the one-step law of the future given the present. Let RX ∼ µ be the present state, X ′ ∼ p(· | X) the coupled future, and µ′ the induced future marginal, dµ′ (x′ ) = p(x′ | x) dµ(x). We assume neither stationarity (µ′ may differ from µ) nor reversibility. The pair (X, X ′ ) has joint law µ̄ on X × X , and this joint law is the sole input — every object below is a functional of µ̄ alone. A coupling is thus any pair of random variables, with the present/future reading as its running instance.

3

The transfer operator. The coupling is carried by its conditional-expectation (transfer) operator T : L2 (X , µ′ ) → L2 (X , µ), (T f )(x) = Ef (X ′ ) | X = x, (1) which sends a future observable to its best present prediction. We also write T = TX,X ′ , when the underlying coupling needs to be specified explicitly. The adjoint of T , T ∗ : L2 (X , µ) → L2 (X , µ′ ) is the backward conditional expectation (T ∗ g)(x′ ) = Eg(X) | X ′ = x′ , and constants are fixed, T 1 = 1 = T ∗ 1. Regularity and the singular system. The single standing assumption in this paper is that T is Hilbert– Schmidt, tr(T ∗ T ) < ∞. This is also the standing assumption of the variational Koopman literature, where it underwrites the singular value decomposition and likewise excludes deterministic dynamics [Wu and Noé, 2020]. A Hilbert–Schmidt operator is compact, so T admits a singular value decomposition (σj , ϕj , ψj )j≥0 [Hsing and Eubank, 2015], T ψj = σj ϕj , T ∗ ϕj = σj ψj , (2) with singular values 1 = σ0 ≥ σ1 ≥ · · · → 0, present functions {ϕj } orthonormal in L2 (µ), future functions {ψj } orthonormal in L2 (µ′ ), and the trivial leading triple σ0 = 1, ϕ0 = ψ0 = 1. The trace is the sum of the P 2 squared singular values, and is the squared Hilbert–Schmidt norm: tr(T ∗ T ) = ∥T ∥HS = j≥0 σj2 . In addition, the trace is related to the χ2 dependence between X and X ′ — the χ2 divergence of the joint law from the independent product of the marginals, written χ2 (X; X ′ ). Let P be the joint law of (X, X ′ ) and r = dP/d(µ ⊗ µ′ ) its density ratio against the productRof the marginals. Then the transfer operator T = TX,X ′ is the integral operator with kernel r, (T f )(x) = r(x, x′ )f (x′ ) dµ′ (x′ ) and we have Z 2 ∥T ∥HS = r2 d(µ ⊗ µ′ ) = 1 + χ2 (X; X ′ ). (3) 2

We thus use ∥T ∥HS and 1 + χ2 (X; X ′ ) interchangeably throughout the paper.

4

Framework and Results

4.1

General Embedding Results

Induced couplings. Let ϕ : X → Y be a measurable feature map. With Y = ϕ(X) and Y ′ = ϕ(X ′ ) the pair (Y, Y ′ ) is again a coupling, the one induced by ϕ, and it carries its own transfer operator TY,Y ′ . The algebra objective. We will be interested in embeddings ϕ that maximize the size of the induced operator, Tϕ := TY,Y ′ . By (3), applied to the induced coupling,
2 ∥Tϕ ∥HS = 1 + χ2 ϕ(X); ϕ(X ′ ) ,

(4)

the predictive χ2 dependence — the total principal inertia — of the embedded present and future. We call (4) the algebra objective because its value depends on ϕ only through the σ-algebra σ(ϕ) the coordinates generate:1 any embedding computing the same information receives the same value, and every product, power, or other measurable function of the coordinates is credited for free. Theorem A.1 in Section A states the relation between (4) and the space Vϕ of Section 1 precisely: Tϕ is the compression of T onto Vϕ and its future counterpart. The class of maps, and the optimum over it. The objective is defined for every measurable ϕ, but what is attainable depends on which maps one is willing to search. We therefore carry the class as an explicit parameter: Φ is a set of measurable maps X → Rm , and ϕ∗ denotes a maximizer of the objective over it, 2

ϕ∗ ∈ arg max ∥Tϕ ∥HS ,

(5)

ϕ∈Φ

1 A function g is σ(ϕ)-measurable iff it is of the form g = ρ ◦ ϕ for some measurable ρ : Y → R — iff it is a function of the

coordinates.

4

2

attaining the value ∥Tϕ∗ ∥HS . Since the criterion sees ϕ only through the algebra it generates, ϕ∗ is determined only up to that algebra, which is the object actually selected. Two properties hold for every Φ: the value 2 never exceeds ∥T ∥HS = 1 + χ2 (X; X ′ ), by the data-processing inequality, and it is non-decreasing as Φ grows, in particular when a coordinate is adjoined to every map in it. Three classes appear below — all measurable maps, the Lipschitz maps that Theorem 4.1 produces, and the networks on the observed coordinates that we actually fit. Remeasuring the state. Suppose the state is observed through a map h : X → Z, so that models must be built from Z = h(X) rather than from X itself. For a general h this loses dependence, by data processing. When h is invertible, however, nothing is lost: the induced coupling (h(X), h(X ′ )) is the original coupling relabeled, and every embedding ϕ of the state translates to the embedding ϕ̃ = ϕ ◦ h−1 of the observation — the same coordinates, now computed from Z, provided the class of maps on Z is rich enough to contain ϕ̃. We return to this in the experiments (Section 5.4). See also Section F. A budget of m coordinates. With Φ a family of maps X → Rm , how large should the dimension m be? We first note that this depends on the complexity of Φ: a single coordinate, m = 1, made complex enough, can be injective, and its generated algebra then carries everything, so the budget m by itself is not an obstacle. The real question is — what m suffices with maps of controlled complexity? For instance, what are we guaranteed for a system with smooth dynamics in a d-dimensional state space and with reasonably behaved maps ϕ? To answer this question, we first define the intrinsic dimension of a system in terms of quantities governed only by the coupling (X, X ′ ) itself. Specifically, we set d∗ to be the box dimension [Falconer, 2014] of the state space X under the predictive diffusion metric of Coifman and Lafon [2006], in which two states are close when their predictive laws are close. This definition parallels classical constructions in dynamics, for instance [Farmer et al., 1983], [Sauer et al., 1991] — with the attractor replaced by the coupling’s predictive geometry, as in Bittracher et al. [2018]. The intrinsic dimension behaves as expected on smooth systems: on a compact d-dimensional state space with smooth enough, nondegenerate transitions we have d∗ ≤ d (Theorem H.2), and the system of Section 5 has d∗ = 10. The full details of the construction are given in Section I. In the next theorem we then show that m = ⌊2d∗ ⌋ + 1 always suffices. This holds no matter how many significant modes the operator carries: the guaranteed budget is set by the dimension of the dynamics, not by the size of its spectrum. Theorem 4.1 (A finite budget suffices). Assume the coupling is Hilbert–Schmidt with finite intrinsic dimension d∗ < ∞. Then for every integer m > 2d∗ there exists a map ϕ : X → Rm , Lipschitz with respect to the predictive metric, whose generated algebra carries the whole operator, 2

2

∥Tϕ ∥HS = ∥T ∥HS = 1 + χ2 (X; X ′ ).

(6)

In particular ⌊2d∗ ⌋ + 1 coordinates suffice. The proof sketch and the full proof are in Section I. The proof builds on general embedology results in Banach spaces [Hunt and Kaloshin, 1999], and is essentially an existence proof, exploiting a certain kind of random linear mappings, and does not otherwise provide a recipe for constructing the map. Maximizing (4) may be viewed as the learning of such a mapping.

4.2

Connection with Other Objectives

Span versus algebra: what is scored. The classical rank-based objective, VAMP [Wu and Noé, 2020] — the objective of VAMPnets [Mardt et al., 2018] — searches for the map ϕ(x) = (ϕ1 (x), . . . , ϕm (x)) ∈ Rm for which T restricted to the linear span of the coordinate functions, span {ϕ1 , . . . , ϕm }, has the largest norm. In Theorem A.1 (Section A) we show that the objective (4) is precisely this criterion with the span replaced by the whole σ-algebra the coordinates generate — products, powers, and every other measurable function of the ϕi included. The span credits only the coordinates themselves — every product and power must be bought with a further coordinate — while the algebra credits them for free. We take the larger algebra, which allows the map ϕ to extract more information at the same budget. 5

Estimation via the variational form. Next, one way to estimate a χ2 divergence is via its variational form. For the coupling (X, X ′ ), with joint law P and independent product of marginals µ ⊗ µ′ as in Section 3, the χ2 is a supremum over functions g ∈ L2 (µ ⊗ µ′ ) of the pair, called critics: Z   Z sup 2 g dP − g 2 d(µ ⊗ µ′ ) = 1 + χ2 (X; X ′ ) (7) g

[Kanamori et al., 2009, Nguyen et al., 2010]. Proofs and attribution are in Section B. Bilinear versus full critics: what is estimated. Applied to the induced coupling (ϕ(X), ϕ(X ′ )), the supremum in (7) computes the χ2 in (4), with critics g(y, y ′ ) now functions of the embedded pair. We use this approach for the experiments in Section 5. Notably, the critic class is where the span–algebra distinction re-enters: restricting the critics g in (7) to a smaller class captures a smaller part of the dependence. Over the full critic class the supremum is the entire 1 + χ2 (ϕ(X); ϕ(X ′ )). Restricting to bilinear critics — functions affine in each argument, g(y, y ′ ) = α + a⊤ y + b⊤ y ′ + y ⊤B y ′ — admits only the coordinates themselves, and the supremum collapses to exactly the VAMP-2 span score of ϕ (see Turri et al., 2026, HaoChen et al., 2021). Finally, we note that the fact that the χ2 captures the full dependence of the coupling is well known. The novelty in this paper is in maximizing ϕ against this criterion rather than the bilinear one, in the motivation for doing so, and in its implications.

5

Experiments

We study masking and the algebra objective on a multicomponent diffusion assembled from published benchmark systems, with a coupling dial κ controlling the interdependence of its components. The system is intrinsically ten-dimensional, and its components are of two kinds: two slow ones with rich spectra, and six fast ones whose modes lie below them. We show that span methods mask its fast components up to large rank — through k = 80, with the transition at k ≈ 100, the product-mode count — while a ten-dimensional maximizer of the algebra objective recovers every component (Section 5.3). Both objectives are invariant under a large class of reparametrizations (Section 4.1), so constructing the representations on the native state is no loss of generality; we test this directly by rerunning the program through a high-dimensional, strongly non-isometric warp of the state (Section 5.4), where an interesting optimization phenomenon also appears. Finally we ask what a representation buys over regressing each target directly from the data: learned from unlabeled observations, the embedding supports prediction from far fewer labeled samples (Section 5.5).

5.1

The Composite System

Components. The system is built from slow and fast components, which play opposite roles in the masking geometry: the slow ones carry rich spectra, so their products supply the leading modes of the joint operator, while each fast one contributes a single mode, placed below those products. The slow part consists of r = 2 independent copies of the lemon-slice diffusion — overdamped Langevin dynamics in the circular multi-well potential of Bittracher et al. [2018], as used by Klus et al. [2020b] and named in Klus et al. [2020a] — p

2 Vcool (x, y) = cos a atan2(y, x) + 10 x2 + y 2 − 1 , with a = 10 wells and diffusion coefficient D = 0.25. The fast part is a six-dimensional hypercube of independent double-well bits, Vhot (x) = 2(x2 − 1)2 per coordinate at D = 0.75, the benchmark family of iVAMPnets [Mardt et al., 2022]. We call the ring components cools and the bits hots, and the state is 2 × 2 + 6 × 1 = 10 dimensional. All experiments use transition pairs at lag τ = 2. Coupling. A dial κ ≥ 0 couples the components: a Kuramoto torque ±κ sin(θ1 − θ0 ) between the two cool component angles, and a directed tilt −κ cos(θ0 − ψj ) xj on bit j, with phases ψj = 2πj/6, so cool0 modulates every bit while no bit feeds back (Figure 1). At κ = 0 the components are independent and the cool marginals are exact nearest-neighbor ring walks; we use κ ∈ {0, 0.5, 1, 2}, chosen so that every component remains live and autonomous at every setting. 6

cool0 (ring, a = 10)

cool1 (ring, a = 10)

±κ sin(θ1 − θ0 )

θ0

θ1

−κ cos(θ0 − ψj ) xj

x0

x1

x2

x3

x4

x5

Figure 1: The coupling at dial κ: a symmetric torque between the two rings, and a phase-staggered tilt from cool0 on each bit, with no feedback. At κ = 0 every arrow vanishes. 1

hot block (6 bits)

σ

∘ MSM (measured) cool marginal (ring walk)

0 1

2

3

4

5

6

7

8

9

σ

1.0 0.8 0.6

cool products

1

25

50

hot-bearing

75

100

100+

joint spectrum rank (closed form, from the marginals)

Figure 2: The composite system’s spectrum. Top: the marginal of one cool, a ring walk, with the measured Markov state model overlaid, and the six bits as a degenerate block at σ ≈ 0.55. Bottom: the joint spectrum at κ = 0, the product of the marginals. The rank axis is broken because the crossover is near rank 100 but not pinned there, the closed form idealizing each cool as a ten-state ring walk. Spectrum, and why masking is predicted. The leading singular values of a single cool component can be approximated by those of a cyclic walk on its wells [Sarich et al., 2010, Prinz et al., 2011], λi = 1 − p(1 − cos(2πi/a)) at hop probability p = 0.080 [Levin and Peres, 2017], yielding values 0.985/0.945/0.895 for the leading nontrivial modes, verified in addition by a Markov state model approximation on simulated data (Figure 2). Each bit contributes a single mode at σ ≈ 0.55. At κ = 0 the joint singular values are products across components, so the whole joint spectrum follows from these marginals by multiplication, with no estimation involved. Figure 2 plots it: the two cools alone generate ar = 100 product modes, and the entire hot block falls below them, the lowest cool product sitting at (1 − 2p)2 ≈ 0.71 against σ1hot = 0.561 for the highest mode that involves a bit. The joint rank saturates at 100. A rank-k span optimum therefore fills with cool products and excludes every bit until k reaches the product count — roughly one hundred coordinates for a system whose intrinsic dimension is ten — while six individually predictable components sit outside the span.

5.2

Estimators and Protocol

Span. VAMPnets [Mardt et al., 2018, Wu and Noé, 2020] train a k-output network on the VAMP score. We train to convergence under a fixed protocol: two restarts, checkpoints to 2400 epochs, restart selection and snapshotting on held-out VAMP-E. Convergence is not a formality here: at every setting we tested, an 7

bit recovery (R2 ) κ=0 κ = 0.5 κ=1 κ=2 floor (cools only) 0.00 0.02 0.12 0.39 span, k = 30 0.00 0.02 0.12 0.39 bilinear critic 0.04 — 0.12 — algebra, m = 10 0.24–0.27 0.24–0.28 0.26–0.32 0.47 ceiling (raw state) 0.26 0.27 0.31 0.48

Table 1: Bit recovery along the coupling dial, mean over the six fast components. The floor is what the cools alone reveal about them. The bilinear-critic row is our own objective at m = 10 with the critic restricted to the span, run at two couplings. Algebra entries are min–max over three, two, three and two runs. under-trained VAMPnet carries transient bit signal that its own objective removes as the held-out score improves (Section 5.4), so a shorter protocol reports spurious recovery. An undersized network produces the same artifact: unable to build the leading product modes, it spends those coordinates on the bits instead (Section D.2). Algebra. We maximize (7) by joint ascent over an m-output embedding ϕ and an unrestricted MLP critic ˆ Restricting the same functional to bilinear g, with snapshot early stopping on the held-out objective value J. critics recovers the span case (Section 4.2). Run in our own machinery, that restriction is indistinguishable from VAMPnet on every readout, which confirms the phenomenon through an independent implementation (Table 1). Representations train on n = 30,000 unlabeled pairs. Writing wℓ for a fully connected trunk of ℓ hidden layers of w units each, the algebra trunk is 642 on the native observation and 1282 on the warped one — never wider than the span trunk it is compared against, which we sweep to 5122 . Readout. Every method is scored the same way: uniform prediction targets — each component’s future state coordinates, identically for every component — regressed from the learned features by one neural readout on separate splits, reported as per-component test R2 . Two reference rows apply that same readout to inputs that are not learned. The first is the full state. Every representation is a function of the state, so under a common readout none can predict a target better than the state itself does: this row is the ceiling, the best score any method could reach. The second withholds the bit coordinates and reads out from the cools alone. This is the floor benchmark : once the components are coupled the cools carry some information about the bits, so a method clears that level without representing a bit at all. It is zero at κ = 0 and rises with the coupling.

5.3

Masking and Recovery

Figure 3 shows bit recovery against span rank k, at the largest trunks we train (5122 ), one curve per κ, with each curve’s floor benchmark and ceiling marked, and the algebra objective at m = 10 as a level marker. The pattern is the same at every coupling: the span sits exactly on its floor benchmark through k = 80 — at κ = 0 that floor is zero and the converged span returns −0.01 at k = 80 across two restarts, at κ = 1 VAMPnet reproduces the floor value 0.117 to the third decimal at both k = 50 and k = 80 — then admits the bits gradually across k = 100–120, near the product count. We call that rank the threshold : the smallest budget at which a span carries anything of the bits’ own content. It does not move with κ: coupling raises the floor (cool-borne bit information grows until, at κ = 2, the cools alone reach 82% of the bits’ predictable content) but not the rank at which the bits’ own content enters. Table 1 gives the same comparison numerically, at one rank per coupling, where the floor and the ceiling are easier to read off than from the figure: the span reproduces the floor, and the algebra objective at m = 10 reaches the ceiling, at every coupling. Cool recovery is uniformly high (0.91–0.96) for every method at every k and κ.

5.4

Warped Observations

Native dynamical coordinates are rarely observed in practice: what is given instead is a larger set of correlated measurements that describe them, which is why a dimension reduction is needed at all. We model this situation by observing the state through a fixed random injective warp S : R10 → R100 — an orthonormal 8

bit recovery (R 2, mean)

laundering floors algebra m = 10 (per κ)

0.5 0.4

κ=2

ceilings

0.3 0.2

κ=1

0.1 κ = 0.5 κ=0 0.0 10

30

50

80

100

120

span rank k

Figure 3: Masking along the coupling dial. Bit recovery of the rank-k span optimum stays on the floor benchmark until k approaches the product-mode count (≈ 100), at every coupling strength, while the algebra objective recovers the bits at m = 10. The floor rises with κ as the cools come to carry more about the bits; the threshold does not move. injection followed by four random affine coupling layers of a normalizing flow — smooth, exactly invertible, and strongly non-isometric (local scale distortion ≈ 7 measured on data pairs). Thus, in this subsection and the next, every method receives the 100-dimensional warped observation as input, in place of the ten-dimensional native state used above. Since the warp is an invertible remeasurement of the state, the achievable values of both objectives are unchanged (Remeasuring the state, Section 4.1), so the population picture is unmoved. The measurements agree: VAMPnet stays masked through k = 50 and admits the bits across k = 80–120, as it does natively. What the warp costs is estimation: past the threshold the best held-out span scores run about a third below native. An optimization phenomenon, and its repair. Optimizing the algebra objective on the warped data by direct joint ascent over ϕ and g falls short: it plateaus at Jˆ ≈ 93 against ≈ 288 on the native state, recovering the bits only partially at full readout (0.19–0.23 against a 0.26 ceiling) and weakly in the few-label regime. The shortfall is not capacity — enlarging ϕ in width or depth strictly lowers the attained value, and enlarging the critic leaves cold features’ scores unchanged — and it is not the estimand, since a diagnostic ϕ frozen at a supervised inverse of S scores at the native level. It is the ascent itself. The repair is standard unsupervised pretraining [Erhan et al., 2010]: an autoencoder on the unlabeled warped stream, trained on reconstruction alone with no lagged pairs and no dynamical term (reconstruction R2 = 0.85, so an imperfect chart suffices), whose encoder warm-starts the ascent. It sees no dynamics and uses no labels, so it remains an optimization aid rather than a second model. Started this way, the ascent restores native-level recovery (bits 0.253 at full readout, three seeds, negligible seed variance). The failure itself is interesting: a 30k-parameter network with n = 30,000 pairs and a smooth invertible target is a regime in which direct optimization would be expected to succeed, and it does not. Our measurements place the difficulty in the joint ascent rather than in the objective or the architecture, and we think identifying it is a worthwhile question in its own right. The objective, not the optimizer, removes the bits. The same warm start separates the two objectives cleanly, and it is a neutral starting point precisely because the chart is dynamics-blind: having never seen a lagged pair, it carries no preference among components into the comparison, so every divergence below is attributable to the objective that follows it. Given the identical autoencoder chart, the free-critic ascent keeps the bits (0.253); the bilinear-critic ascent removes them (−0.01); and a VAMPnet whose trunk is initialized from the chart starts with readout-verified bit information (R2 = 0.244 at its first checkpoint, at the recovery level of the algebra ϕ) and trains it away while its held-out score improves monotonically, reaching −0.01 by convergence (Figure S2). Every convergence trace we recorded, in every setting, shows this decay of early

9

from 300 labels cools bits direct regression 0.40–0.44 0.03 ± 0.04 span, any rank 12–120 — ≤0 autoencoder chart alone 0.49 0.09 ± 0.02 algebra, m = 10 0.77–0.80 0.13 ± 0.02 ceiling (full sample) — 0.26

Table 2: Prediction from few labels on the warped observation, mean ± sd over twenty paired label draws that resample the labels only, with nothing retrained. Spans are at or below zero on the bits at every rank we measured, k = 12 to 120, including the ranks whose full-sample readout does recover them. transient bit signal; the warm-started run certifies that what decays is real, usable information, removed because the span optimum excludes it. This is also why the convergence protocol of Section 5.2 is load-bearing: any under-trained span representation exhibits transient “recovery” that vanishes at convergence.

5.5

Prediction from Few Labels

Why learn a representation at all, when any target can be regressed directly from the raw data? We show that on the warped observation the algebra objective learns representations that support few-shot learning. Representations are trained on n = 30,000 unlabeled pairs, each prediction target is then learned from 300 noisy labeled samples, and we report test R2 against clean targets, resampled over twenty independent label draws with all methods paired on identical draws. Five representations are compared under that protocol: the warped observation itself, a rank-k span, the autoencoder chart, the algebra objective at m = 10, and the native state, which is the ceiling. Table 2 gives the comparison. The algebra representation beats direct regression on both kinds of component, and beats every span rank on the bits. Its advantage is also not inherited from the unsupervised warm start that initializes the ascent: the autoencoder chart alone already improves on direct regression, and training the algebra objective from that chart improves on it again. Span representations sit at or below zero on the bits at every rank k = 12–120, in every setting and at every capacity we measured — including ranks past the threshold whose full-sample readout does recover the bits (native 5122 , k = 120: 0.17–0.19 at full sample, −0.05 at 300 labels even with clean labels). Admission into the span at full sample does not transfer to the few-label regime: the information is spread across ≈ 100 coordinates, and 300 labels cannot reassemble it, while the algebra objective concentrates it in ten. The scope of the advantage is two-fold. Dominant structure is recoverable from few labels through either compression at sufficient rank or the algebra objective. The weak components are recoverable only through the algebra representation — at full sample by any route that does not compress them away, at few labels by nothing else we measured.

6

Conclusion, Limitations, and Future Work

We proposed a criterion for reducing a dynamical system to few coordinates: score the algebra the feature map generates rather than the span of its coordinates, so that a component costs its generators and not its modes. The objective is a χ2 dependence between the embedded present and future, optimized directly by standard density-ratio estimators, and its variational form places VAMP and our criterion in one family, separated only by the critic class — bilinear for the span, unrestricted for the algebra. Experimentally, span objectives exhibit the motivating failure, linear masking — whole components dropped at every rank below the count of interaction modes, a count exponential in the number of components — while the algebra objective recovers them at one coordinate per degree of freedom, and its representations support prediction from few labels where regression on the observation fails. These results extend the operator-theoretic program rather than displace it: what changes is the price of a component — degrees of freedom rather than modes. The warped-observation study (Section 5.4) also yields a finding: from an unsupervised warm start — standard pretraining, blind to the dynamics — the ascent reaches its native value, and the measurements localize the cold-start gap to the joint optimization rather than to the objective. We leave this as future work. Finally, the experimental construction we use is a natural extension of studied benchmark systems, but

10

the data are simulated rather than real — the next step is real data, where the components are not known in advance.

References Galen Andrew, Raman Arora, Jeff Bilmes, and Karen Livescu. Deep canonical correlation analysis. In International Conference on Machine Learning (ICML), pages 1247–1255, 2013. Francis R. Bach and Michael I. Jordan. Kernel independent component analysis. Journal of Machine Learning Research, 3:1–48, 2002. Andreas Bittracher, Péter Koltai, Stefan Klus, Ralf Banisch, Michael Dellnitz, and Christof Schütte. Transition manifolds of complex metastable systems: Theory and data-driven computation of effective dynamics. Journal of Nonlinear Science, 28(2):471–512, 2018. Leo Breiman and Jerome H. Friedman. Estimating optimal transformations for multiple regression and correlation. Journal of the American Statistical Association, 80(391):580–598, 1985. Nicolas Brodu and James P. Crutchfield. Discovering causal structure with reproducing-kernel Hilbert space ϵ-machines. arXiv preprint arXiv:2011.14821, 2020. Steven L. Brunton, Marko Budišić, Eurika Kaiser, and J. Nathan Kutz. Modern Koopman theory for dynamical systems. SIAM Review, 64(2):229–340, 2022. Andreas Buja. Remarks on functional canonical variates, alternating least squares methods and ACE. The Annals of Statistics, 18(3):1032–1069, 1990. Flavio P. Calmon, Ali Makhdoumi, Muriel Médard, Mayank Varia, Mark Christiansen, and Ken R. Duffy. Principal inertia components and applications. IEEE Transactions on Information Theory, 63(8):5011–5038, 2017. Ronald R. Coifman and Stéphane Lafon. Diffusion maps. Applied and Computational Harmonic Analysis, 21 (1):5–30, 2006. Zhijie Deng, Jiaxin Shi, and Jun Zhu. NeuralEF: Deconstructing kernels by deep neural networks. In International Conference on Machine Learning (ICML), pages 4976–4992. PMLR, 2022. Carmeline J. Dsilva, Ronen Talmon, Ronald R. Coifman, and Ioannis G. Kevrekidis. Parsimonious representation of nonlinear dynamical systems through manifold learning: A chemotaxis case study. Applied and Computational Harmonic Analysis, 44(3):759–773, 2018. Dumitru Erhan, Yoshua Bengio, Aaron Courville, Pierre-Antoine Manzagol, Pascal Vincent, and Samy Bengio. Why does unsupervised pre-training help deep learning? Journal of Machine Learning Research, 11:625–660, 2010. Kenneth Falconer. Fractal Geometry: Mathematical Foundations and Applications. John Wiley & Sons, 3rd edition, 2014. J. Doyne Farmer, Edward Ott, and James A. Yorke. The dimension of chaotic attractors. Physica D: Nonlinear Phenomena, 7(1–3):153–180, 1983. Marco Federici, Patrick Forré, Ryota Tomioka, and Bastiaan S. Veeling. Latent representation and simulation of markov processes via time-lagged information bottleneck. In International Conference on Learning Representations (ICLR), 2024. arXiv:2309.07200. Kenji Fukumizu, Francis R. Bach, and Arthur Gretton. Statistical consistency of kernel canonical correlation analysis. Journal of Machine Learning Research, 8:361–383, 2007. Michael J. Greenacre. Theory and Applications of Correspondence Analysis. Academic Press, London, 1984. 11

Arthur Gretton, Olivier Bousquet, Alex Smola, and Bernhard Schölkopf. Measuring statistical dependence with Hilbert–Schmidt norms. In Algorithmic Learning Theory (ALT), pages 63–77. Springer, 2005. Jeff Z. HaoChen, Colin Wei, Adrien Gaidon, and Tengyu Ma. Provable guarantees for self-supervised deep learning with spectral contrastive loss. In Advances in Neural Information Processing Systems (NeurIPS), pages 5000–5011, 2021. Tim Hempel, Mauricio J. del Razo, Christopher T. Lee, Bryn C. Taylor, Rommie E. Amaro, and Frank Noé. Independent Markov decomposition: Toward modeling kinetics of biomolecular complexes. Proceedings of the National Academy of Sciences (PNAS), 118(31):e2105230118, 2021. doi: 10.1073/pnas.2105230118. Tailen Hsing and Randall Eubank. Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators. Wiley Series in Probability and Statistics. John Wiley & Sons, 2015. Brian R. Hunt and Vadim Yu. Kaloshin. Regularity of embeddings of infinite-dimensional fractal sets into finite-dimensional spaces. Nonlinearity, 12(5):1263–1275, 1999. Sungwoo Jeong, Heejun Ryu, Se-Young Yun, and Gregory W. Wornell. Efficient parametric SVD of Koopman operator for stochastic dynamical systems. In Advances in Neural Information Processing Systems (NeurIPS), 2025. arXiv:2507.07222. Takafumi Kanamori, Shohei Hido, and Masashi Sugiyama. A least-squares approach to direct importance estimation. Journal of Machine Learning Research, 10:1391–1445, 2009. Alexander S. Kechris. Classical Descriptive Set Theory, volume 156 of Graduate Texts in Mathematics. Springer, 1995. Stefan Klus, Feliks Nüske, Sebastian Peitz, Jan-Hendrik Niemann, Cecilia Clementi, and Christof Schütte. Data-driven approximation of the Koopman generator: Model reduction, system identification, and control. Physica D: Nonlinear Phenomena, 406:132416, 2020a. Stefan Klus, Ingmar Schuster, and Krikamol Muandet. Eigendecompositions of transfer operators in reproducing kernel Hilbert spaces. Journal of Nonlinear Science, 30(1):283–315, 2020b. Vladimir R. Kostic, Pietro Novelli, Andreas Maurer, Carlo Ciliberto, Lorenzo Rosasco, and Massimiliano Pontil. Learning dynamical systems via Koopman operator regression in reproducing kernel Hilbert spaces. In Advances in Neural Information Processing Systems (NeurIPS), 2022. Vladimir R. Kostic, Karim Lounici, Grégoire Pacreau, Pietro Novelli, Giacomo Turri, and Massimiliano Pontil. Neural conditional probability for uncertainty quantification. In Advances in Neural Information Processing Systems (NeurIPS), 2024a. arXiv:2407.01171. Vladimir R. Kostic, Pietro Novelli, Riccardo Grazzi, Karim Lounici, and Massimiliano Pontil. Learning invariant representations of time-homogeneous stochastic dynamical systems. In International Conference on Learning Representations (ICLR), 2024b. arXiv:2307.09912. H. O. Lancaster. The structure of bivariate distributions. The Annals of Mathematical Statistics, 29(3): 719–736, 1958. David A. Levin and Yuval Peres. Markov Chains and Mixing Times. American Mathematical Society, 2nd edition, 2017. Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton. Deep learning for universal linear embeddings of nonlinear dynamics. Nature Communications, 9(1):4950, 2018. Anuran Makur and Lizhong Zheng. Polynomial spectral decomposition of conditional expectation operators. In 54th Annual Allerton Conference on Communication, Control, and Computing, pages 633–640. IEEE, 2016.

12

Andreas Mardt, Luca Pasquali, Hao Wu, and Frank Noé. VAMPnets for deep learning of molecular kinetics. Nature Communications, 9(1):5, 2018. Andreas Mardt, Tim Hempel, Cecilia Clementi, and Frank Noé. Deep learning to decompose macromolecules into independent Markovian domains. Nature Communications, 13(1):7101, 2022. doi: 10.1038/s41467-02234603-z. Sean Meyn and Richard L. Tweedie. Markov Chains and Stochastic Stability. Cambridge University Press, 2 edition, 2009. Tomer Michaeli, Weiran Wang, and Karen Livescu. Nonparametric canonical correlation analysis. In International Conference on Machine Learning (ICML), volume 48 of PMLR, pages 1967–1976, 2016. XuanLong Nguyen, Martin J. Wainwright, and Michael I. Jordan. Estimating divergence functionals and the likelihood ratio by convex risk minimization. IEEE Transactions on Information Theory, 56(11):5847–5861, 2010. Frank Noé and Feliks Nüske. A variational approach to modeling slow processes in stochastic kinetics. Multiscale Modeling & Simulation, 11(2):635–655, 2013. Sebastian Nowozin, Botond Cseke, and Ryota Tomioka. f-GAN: Training generative neural samplers using variational divergence minimization. In Advances in Neural Information Processing Systems (NeurIPS), pages 271–279, 2016. Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. PyTorch: An imperative style, high-performance deep learning library. In Advances in Neural Information Processing Systems (NeurIPS), 2019. David Pfau, Stig Petersen, Ashish Agarwal, David G. T. Barrett, and Kimberly L. Stachenfeld. Spectral inference networks: Unifying deep and spectral learning. In International Conference on Learning Representations (ICLR), 2019. arXiv:1806.02215. Jan-Hendrik Prinz, Hao Wu, Marco Sarich, Bettina Keller, Martin Senne, Martin Held, John D. Chodera, Christof Schütte, and Frank Noé. Markov models of molecular kinetics: Generation and validation. The Journal of Chemical Physics, 134(17):174105, 2011. Alfréd Rényi. On measures of dependence. Acta Mathematica Academiae Scientiarum Hungaricae, 10(3–4): 441–451, 1959. Jaume Riba and Ferran de Cabrera. Regularized estimation of information via high dimensional canonical correlation analysis. 2020. arXiv:2005.02977. Marco Sarich, Frank Noé, and Christof Schütte. On the approximation quality of Markov state models. Multiscale Modeling & Simulation, 8(4):1154–1177, 2010. Tim Sauer, James A. Yorke, and Martin Casdagli. Embedology. Journal of Statistical Physics, 65(3–4): 579–616, 1991. Matthew S. Schmitt, Maciej Koch-Janusz, Michel Fruchart, Daniel S. Seara, Michael Rust, and Vincenzo Vitelli. Information theory for dimensionality reduction in dynamical systems. arXiv preprint arXiv:2312.06608, 2023. Cosma Rohilla Shalizi and James P. Crutchfield. Computational mechanics: Pattern and prediction, structure and simplicity. Journal of Statistical Physics, 104(3):817–879, 2001. Le Song, Jonathan Huang, Alex Smola, and Kenji Fukumizu. Hilbert space embeddings of conditional distributions with applications to dynamical systems. In International Conference on Machine Learning (ICML), pages 961–968, 2009.

13

Masashi Sugiyama, Taiji Suzuki, and Takafumi Kanamori. Density Ratio Estimation in Machine Learning. Cambridge University Press, 2012. Taiji Suzuki and Masashi Sugiyama. Canonical dependency analysis based on squared-loss mutual information. Neural Networks, 34:46–55, 2012. Taiji Suzuki and Masashi Sugiyama. Sufficient dimension reduction via squared-loss mutual information estimation. Neural Computation, 25(3):725–758, 2013. Naoya Takeishi, Yoshinobu Kawahara, and Takehisa Yairi. Learning Koopman invariant subspaces for dynamic mode decomposition. In Advances in Neural Information Processing Systems (NeurIPS), 2017. Ronen Talmon and Ronald R. Coifman. Intrinsic modeling of stochastic dynamical systems using empirical geometry. Applied and Computational Harmonic Analysis, 39(1):138–160, 2015. Giacomo Turri, Luigi Bonati, Kai Zhu, Massimiliano Pontil, and Pietro Novelli. Self-supervised evolution operator learning for high-dimensional dynamical systems. In International Conference on Learning Representations (ICLR), 2026. arXiv:2505.18671. Aaron van den Oord, Yazhe Li, and Oriol Vinyals. Representation learning with contrastive predictive coding. 2018. arXiv:1807.03748. Johannes von Lindheim. On intrinsic dimension estimation and minimal diffusion maps embeddings of point clouds. Master’s thesis, Freie Universität Berlin, 2018. Lichen Wang, Jiaxiang Wu, Shao-Lun Huang, Lizhong Zheng, Xiangxiang Xu, Lin Zhang, and Junzhou Huang. An efficient approach to informative feature extraction from multimodal data. In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), 2019. arXiv:1811.08979. Matthew O. Williams, Ioannis G. Kevrekidis, and Clarence W. Rowley. A data-driven approximation of the Koopman operator: Extending dynamic mode decomposition. Journal of Nonlinear Science, 25(6): 1307–1346, 2015a. Matthew O. Williams, Clarence W. Rowley, and Ioannis G. Kevrekidis. A kernel-based method for data-driven Koopman spectral analysis. Journal of Computational Dynamics, 2(2):247–265, 2015b. Hao Wu and Frank Noé. Variational approach for learning markov processes from time series data. Journal of Nonlinear Science, 30(1):23–66, 2020.

14

Supplementary Material Contents A. The Projection Lemma. The structural identity the algebra objective rests on: the coupling induced by ϕ is the two-sided compression of T onto the algebra ϕ generates. B. The objective as a χ2 divergence. That the algebra objective is a χ2 divergence, its variational form and the least-squares estimator we use, the critic-class reading that puts span and algebra on one axis, the bilinear ablation, and what a held-out value does and does not certify. C. Additional literature notes. Dependence-maximization methods and the information-bottleneck reductions, component identification, repeated eigendirections, and autoencoder models of the dynamics — several of which carry a cost exponential in the number of components, in coordinates or in batch size. D. Additional experiments. D.1 isolates what the unsupervised warm start supplies on its own. D.2 is the capacity study: an underfit span network does not build the interaction modes its own objective calls for, and may leak bit information in their place. E. The predictive metric. The geometry of the singular system, as against the size that B measures — the χ2 distance between predictive laws, the kernel of T T ∗ and its canonical RKHS, and how the feature class imposes a topology of its own. F. Invariance. The algebra objective is unchanged by any bi-measurable remeasurement of the state, and what that does not give a fitted model. G. What a coordinate buys. Where the algebra objective spends a coordinate, asked twice. Against a rank budget on a product system, in closed form: products and powers are free, a new component is what costs. Against a reconstruction objective: that one allocates by variance, this one by predictability. H. The predictive metrics and the intrinsic dimension. The construction Theorem 4.1 rests on — the forward and backward metrics, the joint embedding, Theorem H.1, and the bound d∗ ≤ d for smooth systems (Theorem H.2). I. Proof of Theorem 4.1. The budget theorem. The definition it depends on is Theorem H.1, in H. J. Experimental details and reproducibility. Compute, seeds, metrics, and the hyperparameter ranges searched.

A

The Projection Lemma

Throughout, write H := L2 (X × X , µ̄) for the joint space, and view every operator in Theorem A.1 as acting on H. We identify L2 (X , µ) with VX via the isometry g 7→ g ◦ pr, (g ◦ pr)(x, x′ ) = g(x) — an isometry because the first marginal of µ̄ is µ — and likewise L2 (X , µ′ ) with VX ′ via f 7→ f (x′ ), the second marginal being µ′ . Under these identifications T is the lift f (x′ ) 7→ (T f )(x), a map VX ′ → VX . For a sub-σ-algebra G write VG for its closed subspace of measurable functions in H and PVG for the orthogonal projection onto it. The present algebra Vϕ ⊆ VX is then the functions of ϕ(X), and the future algebra Vϕ′ ⊆ VX ′ the functions of ϕ(X ′ ). Lemma A.1 (Induced projection). As operators on H, the transfer operator factors into projections T = PVX PVX ′ , and the operator of the coupling induced by any measurable ϕ is the two-sided compression of T — the present projected onto Vϕ , the future onto Vϕ′ , TY,Y ′ = PVϕ T PVϕ′ .

(8)

The proof rests on two elementary facts. (F1) Conditional expectation is orthogonal projection. For a sub-σ-algebra G of the product σ-algebra, PVG h = Eh | G for every h ∈ H: the conditional expectation is the G-measurable function minimizing ∥h − g∥L2 over G-measurable g, which is exactly the orthogonal projection of h onto the closed subspace VG . In particular PVX h depends only on x and PVX ′ h only on x′ . 15

(F2) Nested projections collapse. If A ⊆ B are closed subspaces with orthogonal projections PA , PB , then PA PB = PB PA = PA . Proof of Theorem A.1. First claim, T = PVX PVX ′ . Take h ∈ H. By (F1), PVX ′ h is a function of x′ alone, f (x′ ) := Eh | X ′ = x′ , so PVX ′ h = f ∈ VX ′ , and applying (F1) again, PVX PVX ′ h = PVX f = Ef (X ′ ) | X = x = (T f )(x). Thus PVX PVX ′ restricts to T on VX ′ and annihilates VX⊥′ (where PVX ′ = 0), so the two operators coincide on H. Second claim, (8). The induced pair (Y, Y ′ ) = (ϕ(X), ϕ(X ′ )) is itself a coupling, so the first claim applied to it gives TY,Y ′ = PVϕ PVϕ′ . Because σ(Y ) ⊆ σ(X) we have Vϕ ⊆ VX , and symmetrically Vϕ′ ⊆ VX ′ . Applying (F2) to each nested pair — PVϕ PVX = PVϕ and PVX ′ PVϕ′ = PVϕ′ — with the first claim, PVϕ T PVϕ′ = PVϕ PVX PVX ′ PVϕ′ = (PVϕ PVX )(PVX ′ PVϕ′ ) = PVϕ PVϕ′ = TY,Y ′ , which is (8). Remark A.2. The argument uses only nesting, so the factorization T = PV PV ′ underlying (8) holds for any pair of closed subspaces in place of Vϕ ⊆ VX and Vϕ′ ⊆ VX ′ — in particular for a reproducing-kernel subspace of the algebra, which is what later licenses modeling the operator on an RKHS rather than on all of L2 .

B

The Objective as a χ2 Divergence, and its Direct Estimation

This section carries the proofs and the attribution for Section 4.2: the criterion is a χ2 divergence, its optimal witness is the density ratio r of (9), and least-squares density-ratio estimation computes it variationally, with no covariance inversion and no shared ridge. The criterion is a functional of one object, the density ratio of the coupling. Intrinsically it is the Radon–Nikodym derivative of the joint law P of (X, X ′ ) against the product of its marginals, r =

dP , d(µ ⊗ µ′ )

r(x, ·) =

dp(· | x) , dµ′

(9)

which needs no base measure and lies in L2 (µ ⊗ µ′ ) exactly when χ2 (X; X ′ ) < ∞. If µ, µ′ carry densities pX , pX ′ with respect to a common reference dz, then r is represented by r(x, Rx′ ) = p(x, x′ )/(pX (x)pX ′ (x′ )) = p(x′ | x)/pX ′ (x′ ) — the reference cancels in the ratio — and the integrals (·) dz below are taken against that dz, with µ′ (dz) = pX ′ (z) dz. Its diagonal (singular) expansion [Lancaster, 1958, Rényi, 1959] is X r(x, x′ ) = 1 + σj ϕj (x) ψj (x′ ), (10) j≥1

with {ϕj } orthonormal in L2 (µ), {ψj } orthonormal in L2 (µ′ ), and (σj , ϕj , ψj ) the singular system of T (1). Here σ1 is the maximal correlation and σ0 = 1 on constants. The objective. The criterion of (4) is the total χ2 dependence between the embedded present and future, 2 ∥Tϕ ∥HS − 1 = χ2 (ϕ(X); ϕ(X ′ )), and at ϕ = id it is ZZ X χ2 (X; X ′ ) = (r − 1)2 pX (x) pX ′ (x′ ) dx dx′ = σj2 . (11) j≥1

This scalar is twice the squared-loss mutual information, which is conventionally defined as 12 χ2 , and is the Pearson mean-square contingency and the total principal inertia [Calmon et al., 2017]. It measures the total mass of the singular spectrum. 16

Definition B.1 (χ2 divergence). For probability measures P ≪ Q on a common space, Z  Z  2 dP dP 2 2 χ (P ∥ Q) = − 1 dQ = dQ − 1, dQ dQ

(12)

and χ2 (P ∥ Q) = ∞ if P ̸≪ Q or the integral diverges. It is the f -divergence of f (t) = (t − 1)2 : nonnegative, zero exactly at P = Q, and satisfying the data-processing inequality [Sugiyama et al., 2012]. Lemma B.2 (the objective is a χ2 mutual information). Fix an embedding ϕ, let Pϕ be the joint law of the pair (ϕ(X), ϕ(X ′ )) and Qϕ the product of its marginals. Then
2 ∥Tϕ ∥HS = 1 + χ2 Pϕ ∥ Qϕ ,

(13)

the χ2 analogue of a mutual information: the divergence of the coupling from the independent coupling with the same marginals. Proof. When Pϕ ≪ Qϕ the ratio is rϕ , the object (9) of the embedded pair, and (12) evaluates through the R P 2 expansion (10) to rϕ2 dQϕ − 1 = j≥1 σj (ϕ)2 , which together with the constant mode is ∥Tϕ ∥HS , as in (11). When Pϕ ̸≪ Qϕ the operator is not Hilbert–Schmidt and both sides are infinite. The role of the ratio. By (12) the objective is the squared L2 (Qϕ ) distance of rϕ from the constant 1 — dependence measured as distance from independence — and by (9) the rows rϕ (y, ·) are exactly the predictive laws of the embedded system. So the divergence and the predictor are two functionals of the one object rϕ , and any estimator that recovers rϕ recovers both at once. Least-squares density-ratio estimation. LSIF [Kanamori et al., 2009, Sugiyama et al., 2012] fits r by least squares in L2 (Q) without forming either density. A critic is a function g(x, x′ ) of the pair — it stands in for r, which lives on the product space; its slice R R g(x, ·) at the optimum is the predictive density weight of (9). For any critic g ∈ L2 (Q), using rg dQ = g dP — the defining property of the ratio, which removes r from the criterion, Z Z Z 2 J(g) = 2 g dP − g 2 dQ = r2 dQ − ∥g − r∥L2 (Q) , (14) so sup J(g) = 1 + χ2 (P ∥ Q),

attained at g = r.

(15)

g

This is (13) in variational form — the χ2 case of the convex-dual representation of f -divergences [Nguyen et al., 2010]; the same functional under the name squared-loss mutual information is developed by Suzuki and Sugiyama [2013]. Empirically the first integral runs over the observed pairs and the second over cross-pairings of present with future samples, which sample Q for free. The two consequences used in the body: first, there is no whitening inverse. The critic’s capacity is one budget, but it is spent by the least-squares fit of r itself — concentrated where the dependence is — rather than assigned by the marginal spectrum of the features as the shared whitening ridge of the kernel realization assigns it. (Per-observable regularization is a property of the downstream prediction stage, where each regressed measurable carries its own ridge; the critic does not provide it and is not meant to.) Second, by (14) J(g) ≤ 1 + χ2 for every g, so on a held-out split the estimate approaches the truth from below and cannot be inflated by overfitting. It is not monotone in training — an overfit critic drifts from r and J falls — so the stopping point is itself selected on the held-out value. This is a property of the estimate, not of a critic class: the bilinear (span) critic, whose training is monotone and stable on the decoupled system, develops the same post-peak decline once the components are coupled. The reported model is the held-out snapshot for every critic class alike.

17

The critic-class ablation. Holding the functional, the data, the m = 10 embedding architecture, the optimizer, the schedule and the stopping rule fixed, and changing only the critic class from the free MLP to the bilinear g = 1 + ϕ(x)⊤P ϕ(x′ ), reproduces the span behaviour inside our own machinery. Per-component bit R2 under the neural readout, against a raw-state ceiling of 0.26 at κ = 0 and 0.31 at κ = 1. At κ = 0 the bilinear critic gives 0.037 (cools 0.93), where the free critic gives 0.271/0.242/0.261 over three seeds. At κ = 1 the bilinear critic gives 0.119/0.111/0.119, which is the floor benchmark 0.117 to within seed noise, where the free critic reaches 0.264–0.318. At both couplings the bilinear row sits at the floor to within a few hundredths, far below the ceiling the free critic reaches, so the phenomenon does not depend on the span estimator being a VAMPnet — it follows from the critic class alone. Estimating the dependence. At a fixed ϕ, χ2 (ϕ(X); ϕ(X ′ )) is a divergence of the joint law of the kdimensional pair (ϕ(X), ϕ(X ′ )), an off-the-shelf density-ratio problem. Direct estimators include least-squares importance fitting [Kanamori et al., 2009], squared-loss mutual information and its dimension-reduction use [Suzuki and Sugiyama, 2012, 2013], convex variational divergence estimation [Nguyen et al., 2010, Nowozin et al., 2016], and kernel dependence measures [Gretton et al., 2005], with Sugiyama et al. [2012] the reference treatment. None of these diagonalizes T , so the objective is spectrum-free to evaluate — the singular system enters only the analysis of Theorem 4.1. Remark B.3 (span and algebra are critic classes). Restricting the critic decides which functional (15) computes. Write u = ϕ(X) and v = ϕ(X ′ ) for the present and future features, whitened and mean-free, and
Ĉ01 = Euv ⊤ for their cross-covariance across the lag. A direct computation then gives maxA J 1 + u⊤Av = 2

1 + Ĉ01 — the span (VAMP-2) score of the coordinates — while the supremum over all of L2 (Qϕ ) is (13), F the algebra. A bilinear critic of rank m estimates the top m modes: mode counting re-enters exactly through the critic’s rank, and declining that restriction is the same choice as scoring the algebra in place of the span. Each half is separately on record. The rank-restricted maximum is the variational principle of Wu and Noé [2020]; the same bilinear restriction of (14) recurs as the spectral contrastive loss [HaoChen et al., 2021], in discrete plug-in form in correspondence analysis [Greenacre, 1984, §4.6] and Riba and de Cabrera [2020], and — closest to the present setting — in Turri et al. [2026], who identify its optimal value with the VAMP-2 score for exactly the present–future ratio. Kostic et al. [2024a] fit low-rank models of r with this loss and state that the full-model optimum is the χ2 divergence. Neither endpoint is new and we claim neither. Reading them as one axis is a matter of vocabulary rather than of result: mode counting enters exactly through the critic’s rank, and declining that restriction is the choice of algebra over span. What the reading is for is that it makes the value a quantity to maximize over the embedding. These derivations take the ratio on the full space, where the unrestricted value 1 + χ2 (X; X ′ ) is a constant of the coupling. Ours is the pushforward ratio, whose value 1 + χ2 (ϕ(X); ϕ(X ′ )) moves with ϕ — which is what makes it a criterion for choosing ϕ at all. What the estimate certifies. By (15), restricting g to a class only lowers the value, so for the fitted 2 ϕ the held-out estimate is at most 1 + χ2 (ϕ(X); ϕ(X ′ )) = ∥Tϕ ∥HS , the criterion’s value at that map. The estimate therefore approaches its target from below and cannot be inflated by overfitting the critic. What it certifies is constructive and one-sided: the reported value is attained by a map we exhibit, so m coordinates suffice for it. Showing m coordinates necessary would need an upper bound over all ϕ, which no lower bound supplies and which we do not claim. Convergence of the estimate with the budget. Figure S1 is a diagnostic for the estimator, not a measurement of any dimension. For each m we train the algebra objective on the composite system and ˆ on a grid m = 1–20 at each coupling. Two things are worth reading off it. First, the record the held-out J, absolute scale is set by the critic architecture — values rise by roughly 10% per doubling of the critic width — so only the shape of each curve is identified, not its height, and no comparison against the closed-form ceiling of Theorem G.1 should be drawn from it. Second, the value achievable at budget m is non-decreasing in m, since adjoining a coordinate can only enlarge Vϕ , so any decline past a plateau is estimation error rather than signal. Growth stabilizes on every curve, and past that point the curves move by no more than the spread between optimizer restarts at a single m, so what remains there is run-to-run scatter rather than a trend.

18

d ∗ = 10

300

frontier J ̂

250 200 150 100 50

κ=0 κ = 0.5

κ=1 κ=2

0 1

5

10

15

20

embedding dimension m

Figure S1: Convergence of the objective with the coordinate budget, on the composite system of Section 5. Held-out Jˆ against embedding dimension m at each coupling, seed spread shaded. The vertical scale is set by the critic architecture and is not comparable across architectures, so the curves identify shape and not height. Because the achievable value is non-decreasing in m, any decline past a plateau measures the estimator’s error rather than a loss of signal. Growth stabilizes on every curve, later at stronger coupling. This is an instrument diagnostic: the flattening reports that our estimator stops gaining, not that the coupling has a dimension. Where growth stabilizes moves with the coupling, later at stronger coupling. The plateau itself falls as κ rises, the criterion being a property of the coupling, and coupling spending part of it.

C

Additional Literature Notes

χ2 / dependence-maximization methods. Maximizing a χ2 dependence to select features is an established program. Squared-loss mutual information drives sufficient dimension reduction [Suzuki and Sugiyama, 2013] over a projection that is linear, static, and aimed at an external target — it recovers the central subspace and carries no dimension notion. Principal inertia components [Calmon et al., 2017] and Soft-HGR [Wang et al., 2019] maximize the same dependence over k feature pairs and recover the top-k modes, and the spectral contrastive loss [HaoChen et al., 2021] characterizes its optimum as the top eigenfunctions — in our vocabulary, span readings, which count components by their modes. Our objective instead credits every function of the coordinates, so once a component’s generators are held, its products and powers are free. Schmitt et al. [2023] share a principal aim with our approach — reducing a system by how predictive its variables are of the future — but measure predictability by a mutual information against a bit-rate budget, through an information bottleneck. That objective is considerably harder to compute: solved exactly it needs the full conditional law p(xt+∆t | xt ), which they note is difficult to estimate in practice, so their variational form replaces it with a prior on the code and a contrastive estimate of the predictive information. A similar objective, likewise built on mutual information, appears in Federici et al. [2024], who maximize the mutual information between the embedded present and the embedded future, and estimate it contrastively as well. Their bottleneck term is not directly computable either, and is replaced by a bound requiring an explicit conditional density model. Both therefore inherit the ceiling of the contrastive estimate. INCE is bounded by log B at batch size B [van den Oord et al., 2018], so representing I nats of dependence requires B ≥ eI , and since the information in multicomponent systems typically grows roughly linearly in the number of components (information is 19

additive for independent components), the batch size required grows exponentially in their number. Our objective carries no such bound: at any batch size the held-out value is an unbiased estimate of the criterion at the fitted critic, and the supremum over critics is 1 + χ2 exactly (7). Component identification. The independent-Markov-decomposition line [Hempel et al., 2021] and its deep successor [Mardt et al., 2022] build Markov state models of biomolecular complexes. As discussed in Section 1, the state complexity of a product system grows with the number of components, and they approach this by decomposing the system into weakly coupled subsystems. Repeated eigendirections. Dsilva et al. [2018] observe that a diffusion-maps embedding returns coordinates parametrizing a direction that an earlier coordinate already parametrizes — repeated eigendirections, their instance being cos 2x against cos x — which overstates the dimension of the data and pushes a genuinely new direction below the repeats in the spectrum. These repeats are exactly what the algebra objective credits for free: once cos x is held, every function of it is available and (4) pays nothing further for cos 2x. In contrast to our approach, their solution is not to seek a better embedding but to filter the one computed: each eigenvector is fitted from its predecessors by a local linear regression, and is discarded if that fit succeeds. Every kept direction must therefore first appear in the computed spectrum, and on the system of Section 5 the first fast component appears only past rank one hundred — a hundred eigenvectors and a hundred such tests to reach it. That count is the product mode count, so it grows exponentially in the number of slow components: the high-rank spectral computation is paid in full and only the final dimension is reduced, where (4) never forms a spectrum at all (Section B). Autoencoder methods. Koopman autoencoders [Lusch et al., 2018, Takeishi et al., 2017] linearize the dynamics with a finite-dimensional linear model. That requirement implies that a rank objective is the optimal one, so here the rank restriction is structurally necessitated. We require less: no model of the dynamics is fit at all, the embedding only generates the algebra that (4) scores, and predictions are made from it afterwards by regression. Finally, we note that we use an autoencoder in Section 5.4, for the purposes of pre-training. This autoencoder does not model the dynamics: it is trained to reduce the dimension of present-only data, with no future and no lagged pairs, and is used as a warm start for our objective. Adding the dynamic part, our objective itself, then improves on it (Section D.1).

D

Additional Experiments

D.1

The Warm Start Is Not the Advantage

The few-label results of Section 5.5 use an algebra ϕ trained from an autoencoder warm start (Section 5.4), which raises the question of how much of the advantage the unsupervised chart already supplies. We answer it by reading out the chart itself, with no algebra training on top: the same architecture, the same checkpoint the algebra runs start from, the same twenty paired label draws and the same readout. At 300 noisy labels on the warped observation, averaged over three seeds: representation direct regression autoencoder chart alone algebra ϕ from that chart

cools 0.42 0.49 0.78

bits 0.030 ± 0.038 0.089 ± 0.019 0.134 ± 0.017

The chart alone carries real information about the fast components — three times what direct regression recovers — so part of the gap to direct regression is due to unsupervised pretraining rather than to the criterion. The remainder is due to the criterion, and it is the larger part on the dominant components: training the algebra objective from that chart adds +0.29 on the cools and +0.045 on the bits, winning on all twenty paired draws in both cases (exact two-sided sign test, p = 1.9 × 10−6 ). The comparison is exact in the sense that the two rows differ only in whether the objective was applied.

20

algebra ϕ, same chart

22

bits R 2

0.2 0.1

20

0.0

18 200

600

1200

held-out VAMP-E

0.3

2400

training epoch

Figure S2: Removal is the optimization. Initialized from the same unsupervised chart, the span objective trains away readout-verified bit information (orange, left axis) as its own held-out score improves (blue, right axis), while the algebra objective from that chart retains it. Two restarts.

D.2

Capacity and the Reachability Wall

An underfit span network does not build the interaction modes its own objective calls for, and leaks bit information into the slots they vacate. Run at a 642 trunk, the sweep of Figure 3 shows the bits entering already at k ≈ 50–80 — a reachability effect, not a sampling one, the span estimator paying twice, once in rank and once in the capacity to compose roughly one hundred interaction functions. That is why Figure 3 is reported at the largest trunks we train. This section records the per-mode evidence. At κ = 0 the product eigenfunctions are available in closed form, so we can measure, per mode, whether the trained network can build them: a linear readout of each product mode from the trained features shows the interaction frontier of a 642 trunk frozen at low order — the cool-cool interaction c3 ×c3 reads out at zero for every k up to 120, while at 2562 –5122 it is built at 60–80% of its ceiling. The sharpest control is a matched pair at essentially identical singular value: the single-component mode c5 ×1 (σ = 0.838) is built by the 642 trunk while the interaction c3 ×c2 (σ = 0.837) is not — the failure is composition, not spectral depth. The same reading comes from the other end of the spectrum: the bits, at σ = 0.55 the lowest modes in play, are built well before interactions that outrank them. The slots the unreachable interactions vacate go to the bits, and the exchange is visible in both directions: down the k-axis at 642 the interaction frontier freezes while bits flood in between k = 50 and 80, and up the width axis at fixed k = 80 the interactions arrive while the bits drain (bit R2 0.25 → −0.01 from 642 to 5122 , which is 0.83 → −0.03 of the ceiling in the fraction-built units Figure S3 plots). Even at 5122 a small unbuildable core remains (the deepest diagonal products), which is why admission in Figure 3 is gradual rather than sharp: the bits enter at the rank the product count predicts, but reach only 0.19 of their 0.26 ceiling by k = 120, a fifth past the wall. The under-capacity and under-training artifacts coincide numerically: at κ = 0, k = 80, the 5122 trunk reads 0.23 at 200 epochs and decays to −0.01 by 2400, the early value coinciding with what a 642 trunk reports at convergence.

E

The Predictive Metric

The objective of Section B measures the size of the singular spectrum. This section concerns its geometry: the predictive metric dS of (23) and its embedding Jf . Like the objective, dS is a functional of the density ratio r of (9) — but unlike the objective it is the singular system. The intrinsic dimension is built from dS but is not its dimension. By Theorem H.1, d∗ is the box dimension of the joint embedding J = (Jf , Jb ) under the joint metric dJ , which separates two states whenever they differ in what they predict or in what they postdict. The forward half alone carries dS and has a dimension of its own, in general a smaller one. Section H assembles dJ from dS and its backward counterpart and defines d∗ there.

21

fraction of mode built

1.2

c5×1 σ.84 c3×c2 σ.84 c3×c3 σ.79 bits σ.55

1.0 0.8 0.6 0.4 0.2 0.0 12

50

80

120

642

1282

2562

5122

trunk width (k = 80)

span rank k (642)

Figure S3: The reachability wall, on both axes. Fraction of each mode built, as a linear readout of that mode from the trained features (R2 /σ 2 ). Left, against span rank at a 642 trunk: the bits are built by k = 80 even though σ = 0.55 makes them the lowest modes plotted, while the interactions above them stay unbuilt. Right, at fixed k = 80 against trunk width: the interactions arrive and the bits drain, the same slots changing hands. c5 ×1 and c3 ×c2 are a matched pair at σ ≈ 0.84, built and unbuilt respectively, so what fails is composition rather than spectral depth. The diffusion metric of (23) is the χ2 distance between one-step predictive laws,
2 Z X
2 p(z | x) − p(z | y) 2 2 dS (x, y) = dz = ∥r(x, ·) − r(y, ·)∥L2 (µ′ ) = σj2 ϕj (x) − ϕj (y) . pX ′ (z)

(16)

j≥1

It admits three readings: the χ2 distance between the predictive laws p(· | x) and p(· | y), the L2 (µ′ ) distance between the rows r(x, ·) of the density ratio, and the weighted Euclidean distance in the singular coordinates Jf (x) = (σj ϕj (x))j . It is the correspondence-analysis χ2 distance between the row profiles of the kernel [Greenacre, 1984] and the diffusion distance [Coifman and Lafon, 2006] of the symmetrized round-trip operator T T ∗ . Unlike the asymmetric χ2 (P ∥ Q), it is symmetric with the fixed reference µ′ , and it vanishes exactly when p(· | x) = p(· | y), that is, when x and y are the same predictive state. The kernel of T T ∗ . Let T ∗ : L2 (µ) → L2 (µ′ ) be the adjoint of T , the backward conditional expectation (T ∗ g)(x′ ) = Eg(X) | X ′ = x′ . The round-trip operator T T ∗ : L2 (µ) → L2 (µ) — step forward x 7→ x′ ∼ p(· | x), then back x′ 7→ x′′ by the posterior p(· | x′ ) — is self-adjoint, positive, and Markov (T T ∗ 1 = 1), hence a reversible Markov operator with stationary law µ and T T ∗ ϕj = σj2 ϕj . Its integral kernel with respect to µ, R defined by (T T ∗ f )(x) = K(x, y)f (y) dµ(y), is Z X p(z | x) p(z | y) K(x, y) = r(x, ·), r(y, ·) L2 (µ′ ) = dz = 1 + σj2 ϕj (x) ϕj (y). (17) pX ′ (z) j≥1

This is the χ2 affinity between the two predictive laws, and its diagonal K(x, x) = 1 + χ2 p(· | x) ∥ µ′ measures how far the future of x departs from the average future µ′ . The predictive metric is exactly the distance induced by this kernel, 2

dS (x, y)2 = K(x, x) − 2K(x, y) + K(y, y) = ∥r(x, ·) − r(y, ·)∥L2 (µ′ ) ,

(18)

so dS is the diffusion distance [Coifman and Lafon, 2006] operator T T ∗ with K as its P of2t the symmetrized 1 2 2 kernel, at diffusion time t = 2 in the convention Dt = j λj (∆ϕj ) with λj = σj2 . Passing to T T ∗ is what 22

lets the embedding machinery for self-adjoint operators apply to the non-reversible coupling T . We avoid the term spectral dimension for d∗ : in the Dirichlet-form literature that name is reserved for the exponent ds = 2α/β read off the on-diagonal heat-kernel decay pt (x, x) ≍ t−α/β , a short-time quantity that a metric taken at a fixed lag cannot see. Remark E.1 (The canonical RKHS of a coupling). Equations (17) and (18) exhibit dS as a reproducingkernel distance. The feature map is x 7→ r(x, ·) = dp(· | x)/dµ′ ∈ L2 (µ′ ), its reproducing kernel is K, and dS (x, y) = ∥K(x, ·) − K(y, ·)∥HK . The space is intrinsic to the coupling: HK = range(T ), the predictable observables f (x) = Eg(X ′ ) | X = x, and ∥f ∥HK is the L2 (µ′ ) norm of the minimal g that predicts f , so less predictable directions (smaller σj ) carry more norm. In particular every predictable observable is dS -Lipschitz: for f = T g, |f (x) − f (y)| = ⟨g, r(x, ·) − r(y, ·)⟩L2 (µ′ ) ≤ ∥g∥L2 (µ′ ) dS (x, y), and the supremum over ∥g∥L2 (µ′ ) ≤ 1 gives the dual form dS (x, y) = sup{ |T g(x) − T g(y)| : ∥g∥L2 (µ′ ) ≤ 1 }, an integral-probability metric [Gretton et al., 2005]. Thus dS is the smallest metric in which T smooths — the coordinates σj ϕj are each 1-Lipschitz and Jf is an isometry onto MT [Coifman and Lafon, 2006] — and this operator-analytic half of Theorem 4.1 is unconditional, leaving all of the reduction’s difficulty in the geometry of MT rather than the action of T . Unlike a method that fixes a kernel on the observation space, nothing is chosen here: the feature r(x, ·) is a Radon–Nikodym derivative, hence invariant under a common bi-measurable reparametrization x 7→ h(x) (Theorem F.1), and K with it. Any characteristic kernel recovers the same canonical system (σj , ϕj , ψj ) after whitening [Fukumizu et al., 2007], so whitening is exactly the passage from a chosen kernel to K, the whitened linear instance being VAMP [Wu and Noé, 2020]. The kernel K itself is classical
(the kernel of T T ∗ ), and it is finite precisely in the Hilbert–Schmidt regime, K(x, x) = 1 + χ2 p(· | x) ∥ µ′ < ∞, degenerating in the deterministic limit where r(x, ·) leaves L2 (µ′ ). Estimating the predictive metric. The metric dS requires the predictive laws p(· | x) or, equivalently, the leading singular functions σj ϕj . These are the estimands of nonlinear canonical analysis: alternating conditional expectations [Breiman and Friedman, 1985], functional canonical variates [Buja, 1990], nonparametric canonical correlation [Michaeli et al., 2016], kernel canonical correlation [Bach and Jordan, 2002, Fukumizu et al., 2007], and conditional mean embeddings or kernel transfer operators [Song et al., 2009, Klus et al., 2020b]. The whitened linear case is the VAMP score [Wu and Noé, 2020], and closed forms are available when the singular functions are orthogonal polynomials [Makur and Zheng, 2016]. The reversible special case recovers the classical diffusion map [Coifman and Lafon, 2006]. The kernel-imported metric. The nearest prior construction of a reduced predictive embedding, the reproducing-kernel ϵ-machine [Shalizi and Crutchfield, 2001, Brodu and Crutchfield, 2020], metrizes predictive laws by a maximum-mean-discrepancy [Gretton et al., 2005] induced by a chosen kernel rather than by dS . That distance depends on the kernel and is not invariant under a bi-measurable reparametrization of the observations, whereas dS is the distance of the canonical kernel K (Theorem E.1) — an arbitrary-kernel proxy in place of the coupling’s own. Because box dimension is a metric — not a topological — invariant, the two give different dimensions, and the intrinsic dimension of Theorem H.1 is the box dimension under the canonical dS . The feature class carries a topology, and it can bias d∗ . The lower bound survives only while every ϕ in the maximization is admissible, that is, continuous in dJ (Theorem H.1). A network on X ⊆ Rd instead enforces continuity in the native metric, and the two classes are in general incomparable. If the singular functions ϕj are native-smooth, a native-continuous ϕ can still separate dS -equivalent states — it is then dS -discontinuous, and the value it attains is not one a dS -continuous map could certify. If instead the ϕj are native-rough or discontinuous, as under a non-smooth observation of the state, the admissible optimum is native-discontinuous and lies outside the network class, which approximates it by several smooth coordinates and over -counts d∗ . Enlarging the class past admissibility therefore deflates the dimension below its invariant value: the operative constraint is fidelity to dS , not network capacity. A shared limitation the algebra softens. This native-continuity prior is common to all neural transferoperator estimators — VAMPnets [Mardt et al., 2018], deep canonical correlation [Andrew et al., 2013],

23

Koopman autoencoders [Lusch et al., 2018] — none of which can represent a native-discontinuous singular function, since the invariance resides in the operator (the σj and the ϕj as L2 objects) and any network reaches it only through the native topology. Relative to these span methods the algebra is the less exposed: a mode method must carry every retained mode as a network output, including the roughest high harmonics, whereas here the network represents only the generators, and their products — the harmonics — are scored through the algebra rather than represented. The network thus carries the fewest and smoothest functions, so the native-smoothness bottleneck binds later than for span methods, and only once the generators themselves become discontinuous.

F

Invariance

Let h : X → X̃ be measurable and injective with measurable inverse on its image, and form the remeasured coupling X̃ = h(X), X̃ ′ = h(X ′ ), with µ̃ = h∗ µ, induced µ̃′ = h∗ µ′ , and operator T̃ = TX̃,X̃ ′ . Pulling functions back along h, Uh : L2 (µ̃) → L2 (µ), (Uh u)(x) = u(h(x)), is an isometry because µ̃ = h∗ µ, and it is onto because h−1 is measurable on h(X ): any v ∈ L2 (µ) is v = Uh (v ◦ h−1 ). Hence Uh is unitary, and Uh′ is its analogue L2 (µ̃′ ) → L2 (µ′ ). Surjectivity of h plays no role, since X̃ \ h(X ) is µ̃-null. When X and X̃ are standard Borel the measurability of h−1 is automatic, injectivity alone sufficing by the Lusin–Suslin theorem [Kechris, 1995]. No surjectivity is required, so observations that embed the state in a larger space, as in the warped and high-dimensional observations of Section 5, are covered. Injectivity is exactly what makes Uh ≤ ∥TY,Y ′ ∥HS , the unitary rather than merely isometric — for a general measurable h one gets T̃Ỹ ,Ỹ ′ HS compression inequality of Theorem A.1. Lemma F.1 (Invariance). For any ϕ̃ : X̃ → Rm put ϕ = ϕ̃ ◦ h, Y = ϕ(X), and Ỹ = ϕ̃(X̃). The induced operators are unitarily conjugate, T̃Ỹ ,Ỹ ′ = Uh−1 TY,Y ′ Uh′ , (19) so T̃Ỹ ,Ỹ ′ = ∥TY,Y ′ ∥HS : the criterion gives ϕ̃ on the relabeled coupling the same value it gives ϕ̃ ◦ h on HS the original. The predictive metric is fixed by the singular system, so the same conjugation carries dJ along. Proof of Theorem F.1. Conditioning on X̃ = x̃ is conditioning on X = h−1 (x̃), so for u ∈ L2 (µ̃′ ), (T̃ u)(x̃) = Eu(h(X ′ )) | X = h−1 (x̃) = (Uh−1 T Uh′ u)(x̃), which is (19) for the full algebras. For the restriction, g is σ(ϕ̃)-measurable iff g = ρ ◦ ϕ̃ for some ρ, and then Uh g = ρ ◦ ϕ is σ(ϕ)-measurable. Conversely every σ(ϕ)-measurable ρ ◦ ϕ equals Uh (ρ ◦ ϕ̃), so Uh Vϕ̃ = Vϕ — this step needs no injectivity, only ϕ = ϕ̃ ◦ h. Since Uh is unitary, by fact (F1) a unitary carries the projection onto a subspace to the projection onto its image, PVϕ̃ = Uh−1 PVϕ Uh and PV ′ = Uh−1 ′ PV ′ Uh′ . Substituting into ϕ ϕ̃

T̃Ỹ ,Ỹ ′ = PVϕ̃ T̃ PV ′ , ϕ̃

−1 T̃Ỹ ,Ỹ ′ = Uh−1 PVϕ Uh · Uh−1 T Uh′ · Uh−1 PVϕ T PVϕ′ Uh′ , ′ PV ′ Uh′ = Uh ϕ

which is (19). Unitary conjugation preserves the Hilbert–Schmidt norm. The predictive metric is a function of the operator’s singular system, which the same equivalence carries along, so dJ -continuity is preserved. Finally ϕ̃ 7→ ϕ̃ ◦ h is a bijection between feature maps on X̃ and on X , with inverse ϕ 7→ ϕ ◦ h−1 — here injectivity is used — so the two problems have the same achievable values at every budget. Remark F.2 (Invariance of the fitted pipeline). Theorem F.1 is a statement about the criterion, not about any estimator. A fitted ϕ is drawn from a class defined by regularity in the observed coordinates, and that class is preserved by h only when h preserves the native topology. Asking h to be a homeomorphism — bi-continuous, not smooth — suffices, and as a bi-measurable map it still leaves dJ and d∗ fixed, so the estimator is equivariant and not only the target. This is also the class on which a kernel fixed on the raw observations fails, ambient smoothness being a metric and not a topological invariant. Section 5.4 tests it.

24

G

What a Coordinate Buys

Recall that the algebra objective scores a map ϕ by the size of the operator of the coupling ϕ induces, 2 ∥Tϕ ∥HS = 1 + χ2 (ϕ(X); ϕ(X ′ )), which is (4). This section evaluates that score exactly in two settings where it can be set against what a different criterion would choose from the same budget. The first asks where a rank budget goes on a system assembled from independent components. The second asks the same of a reconstruction objective. Both pose one question — given one more coordinate, what is it spent on — and the two rivals answer it in the same way as each other. A rank method ranks the candidates by singular value and buys the largest. A reconstruction objective ranks directions by variance and keeps the largest. Each scores a candidate by a magnitude the candidate carries on its own, with no reference to what the coordinates already held can express. The algebra objective (4) scores it instead by what it adds: a coordinate is worth something exactly insofar as it enlarges the algebra σ(ϕ) in a direction carrying present–future dependence. So a mode that is already a product of coordinates held is worth nothing however large its singular value, and a direction that decorrelates within one step is worth nothing however large its variance. The two computations below are those two sentences, made exact.

Products are free, new components are not Q Assume X = (X1 , . . . , Xm ) has independent coordinates and the transition factorizes, p(x′ | x) = i pi (x′i | xi ), so the components evolve without interacting. Let Ti be the operator of component i, with singular values 1 = σi,0 > σi,1 ≥ · · · and left singular functions ϕi,n , and write X 2 2 gi := ∥Ti ∥HS = σi,n (20) n≥0

for its total spectral mass, every order included. Call a component non-degenerate if its leading singular function ϕi,1 generates σ(Xi ), which holds whenever ϕi,1 is injective. Two elementary components qualify. A Gaussian one — Xi ∼ N (0, 1) with (Xi , Xi′ ) jointly Gaussian of correlation ρi — has Hermite singular functions and σi,n = ρin , so gi = 1/(1 − ρ2i ). A two-state one — Xi ∈ {±1} uniform, flipping with probability (1 − λi )/2 — has the single nontrivial mode ϕi,1 = xi at σi,1 = λi , so gi = 1 + λ2i . The second is the idealization of the bits of Section 5. Proposition G.1. Suppose every component is non-degenerate, and let ϕ select a subset I ⊆ [m] of them through their generating coordinates (xi )i∈I . Then Y Y
2 ∥Tϕ ∥HS = gi , i.e. χ2 ϕ(X); ϕ(X ′ ) = gi − 1. (21) i∈I

i∈I

Q 2 Consequently: (1) selecting every component gives the unconditional maximum ∥T ∥HS = i gi ; (2) products are free — adjoining a product ϕi,1 ϕj,1 of two already-selected components leaves the value unchanged, so a coordinate spent on it buys nothing;Q(3) when
σℓ,1 > 0, spending that same coordinate on an unselected component ℓ ∈ / I raises the value by g (gℓ − 1) > 0. i i∈I Proof. The algebra and non-degeneracy gives
objective sees ϕ only through the σ-algebra it generates,
′ ′ σ(ϕ) = σ (Xi )i∈I , so we may take ϕ(x) = (xi )N . The induced pair (Y, Y ) = (X ) , (X ) is then a i∈I i i∈I i i∈I product coupling in its own right: its joint law is i∈I law(Xi , Xi′ ) and its two marginals areQthe corresponding products, so its density ratio against the product of its marginals factorizes, rY (y, y ′ ) = i∈I ri (yi , yi′ ) with ri = dpi (· | xi )/dµ′i . By (3) applied to (Y, Y ′ ) and Fubini, Z YZ Y Y 2 2 2 ′ ∥Tϕ ∥HS = rY d(µY ⊗ µY ) = ri2 d(µi ⊗ µ′i ) = ∥Ti ∥HS = gi . i∈I

i∈I

i∈I

Claim (1) is I = [m]. For (2), ϕi,1 ϕj,1 is a function of (xi , xj ), so adjoining it leaves σ(ϕ), hence the value, unchanged. For (3), adjoining component ℓ replaces I by I ∪ {ℓ} and multiplies the value by gℓ . Only the law of the induced pair enters, so no compression of T and no singular system of the joint operator is needed — the product structure of (Y, Y ′ ) does all the work. The singular values appear only to name gi in (20) and to identify which coordinate a rank method would pick. 25

What the two claims mean. Claims (2) and (3) are one coordinate spent two ways. Suppose ϕ = (x1 , x2 ) and one further coordinate is available. A candidate for it is the product ϕ1,1 ϕ2,1 — not an arbitrary choice, but a genuine singular function of the joint operator, of singular value λ1 λ2 , and when λ1 λ2 > λ3 it is the third-largest mode of the entire system, so a rank-3 method returns exactly x1 , x2 , x1 x2 . Yet that product is already a function of (x1 , x2 ). The generated algebra does not grow, I is still {1, 2}, and the value is still g1 g2 : the coordinate buys nothing whatsoever. Spent on x3 instead it buys g1 g2 (g3 − 1) > 0, strictly positive as soon as the third component has any nontrivial mode. That is the double meaning of products are free. The algebra objective credits every product of the coordinates already held, so a coordinate need never be spent on one — and therefore a coordinate spent on one is spent on something already owned. The smallest instance, and the one we run. Take three two-state components with λ1 λ2 > λ3 . The rank-3 optimum is {x1 , x2 , x1 x2 }, whose algebra is σ(X1 , X2 ) and therefore contains no non-constant function of X3 at all: the third component is not modeled coarsely, it is absent. The algebra objective selects {x1 , x2 , x3 } and exceeds the rank-3 value by (1 + λ21 )(1 + λ22 )λ23 . When λ3 > λ1 λ2 the two selections coincide. So the mode ordering decides whether a rank budget goes wrong, and never affects what (4) prefers. This is linear masking at m = 3, exactly. The system of Section 5 is the same phenomenon with eight components, a budget near one hundred, and measured rather than closed-form quantities — a fuller demonstration, though not a sharper one than this. Remark G.2. A feature violating non-degeneracy — for instance ϕ = x2i , which discards the sign of xi and so captures only the even-order modes of component i — can only lower the value, never exceed the all-component maximum of claim (1).

Autoencoder Reconstruction vs Dynamic Objective An autoencoder model of the dynamics seeks an encoder E : X → Rk , a decoder D : Rk → X , and a latent map M — linear in the Koopman autoencoders of Lusch et al. [2018] — with D ◦ E(x) ≈ x,

D ◦ M ◦ E(x) ≈ x′ .

(22)

The reconstruction requirement makes (22) allocate the budget k by variance rather than by dynamics: minimizing reconstruction error rewards E for keeping the directions of x that are most costly to regenerate, the high-variance ones, whether or not they carry any dynamical information. Write x = x1 ⊕ x2 with x1 high variance but decorrelating in one step — a per-step nuisance — and the predictable dynamics in a lower-variance x2 . Reconstruction spends coordinates on x1 , since dropping it dominates the error, whereas (4) spends none, since x1 contributes no dependence between X and X ′ . In the factorized setting of Theorem G.1 this is exact: a per-step nuisance factor has gi = 1 and zero selection gain in claim (3), at any variance. (A static attribute — x1 frozen, x′1 = x1 — is the opposite case: it is a σ = 1 mode that (4) does retain, being perfectly coupled across time.) Reconstruction methods are in this sense a nonlinear analogue of PCA — a linear autoencoder is PCA — selecting directions by variance, where (4) selects by predictability. This is a property of the loss, not of the decoder, and dropping reconstruction does not repair it. Keeping only D ◦ M ◦ E(x) ≈ x′ — modeling x′ | x — constrains the encoder only on the present support; under a non-stationary coupling µ′ ̸= µ the future occupies a different region, so iterating the model applies E off the support it was fit on. The role of the latent also becomes unclear: a conditional generative model of x′ | x places its latent on the transition noise — the variation in x′ not determined by x — which is not a representation of the state at all. Recovering a representation requires bottlenecking x, which is the one-sided predictive-sufficiency objective once more, now committed to a likelihood we do not need. Finally, even at its optimum the autoencoder recovers only the conditional mean: trained under an L2 loss, D ◦ M ◦ E converges to EX ′ | X = x, a single predicted next state, and says nothing about the spread of x′ around it. The operator T instead encodes Ef (X ′ ) | X for every observable f , and that family determines the full conditional law. The objective (4) carries no decoder and never reconstructs x.

26

H

The Predictive Metrics and the Intrinsic Dimension

This section carries the construction that Theorem 4.1 rests on. The singular system of T puts two metrics on the state space, one forward and one backward. Stacking their embeddings gives the joint embedding whose box dimension is the intrinsic dimension d∗ of Theorem H.1. We then exhibit a class of systems on which d∗ is finite and behaves as expected (Theorem H.2), which is what makes the budget of Theorem 4.1 a finite one. The predictive metrics. The singular system puts a metric on the state space. The forward predictive distance between two states is the χ2 distance2 between the futures they predict,
dS (x, y)2 = χ2 p(· | x) p(· | y) 2

= ∥r(x, ·) − r(y, ·)∥L2 (µ′ ) X
2 = σj2 ϕj (x) − ϕj (y) ,

(23)

j≥1

the weighted Euclidean distance in the left singular coordinates, vanishing exactly when x and y predict the same future. Its time-reversed counterpart, the backward distance, is built identically from the right singular functions and the backward kernel q(· | x) of T ∗ (the previous state given the current),
dS ∗ (x, y)2 = χ2 q(· | x) q(· | y) X
2 (24) = σ 2 ψ (x) − ψ (y) . j

j

j

j≥1

Both are coordinate-free: because r is a Radon–Nikodym derivative, dS and dS ∗ are unchanged when the state is remeasured by any bi-measurable h. The joint embedding and the intrinsic dimension. The forward and backward maps Jf (x) = (σj ϕj (x))j≥1 and Jb (x) = (σj ψj (x))j≥1 embed the state into ℓ2 as the predictive manifolds MT = Jf (X ) and MT ∗ = Jb (X ), carrying dS and dS ∗ [Coifman and Lafon, 2006]. Stacking them gives the joint (bidirectional) embedding J = (Jf , Jb ), an isometry from (X , dJ ) onto J(X ) ⊆ ℓ2 ⊕ ℓ2 under dJ (x, y)2 = dS (x, y)2 + dS ∗ (x, y)2 ,

(25)

which separates two states whenever they differ in what they predict or in what they postdict. Definition H.1 (Intrinsic dimension). The intrinsic dimension of the coupling is the upper box-counting dimension of the joint embedding under the joint metric,
d∗ := dimbox J(X ), dJ , (26) the exponent governing the growth as ε → 0 of the number of dJ -balls of radius ε needed to cover J(X ) [Falconer, 2014]. We take the upper version throughout, which is the one the covering bounds below and the embedding theorem require, and we assume J is defined at every state, that is K(x, x) and K ∗ (x, x) are finite pointwise rather than merely almost everywhere — a standing regularity condition, needed because µ and µ′ may differ and no single null set is exceptional for both. The next proposition exhibits a class of systems with bounded intrinsic dimension: smooth, nondegenerate dynamics on a compact d-dimensional state space. Proposition H.2 (A class of systems with bounded intrinsic dimension). Let (X , ρ) be a compact metric space with upper box dimension d, carrying a finite reference measure dz, and suppose the transition law has a density p(z | x) with respect to dz such that, for constants 0 < c ≤ C < ∞ and L < ∞, 2 The symbol χ2 carries two meanings. Earlier, χ2 (X; X ′ ) in (3) is a dependence between the pair (X, X ′ ), comparing their
joint law with the product of the marginals. Here χ2 p(· | x) ∥ p(· | y) is a divergence between two fixed distributions, the ′ predictive laws of x and y. It is taken against the reference µ rather than the second argument, which makes it symmetric — the χ2 distance of correspondence analysis [Greenacre, 1984].

27

(i) c ≤ p(z | x) ≤ C for all x, z, and (ii) |p(z | x) − p(z | y)| ≤ L ρ(x, y) and |p(x | w) − p(y | w)| ≤ L ρ(x, y) for all z, w and x, y. Then for every initial law µ on X the coupling is Hilbert–Schmidt and d∗ ≤ d. In particular, Theorem 4.1 guarantees that ⌊2d⌋ + 1 coordinates always suffice for such a system. R R Proof. Write |X | = dz for the total reference mass. The future marginal has density pX ′ (z) = p(z | 2 ∈ [c, C], so the density ratio r(x, z) = p(z | x)/pX ′ (z) is bounded by C/c. Hence ∥T ∥HS = w) R dµ(w) 2 ′ 2 ∗ r d(µ ⊗ µ ) ≤ (C/c) < ∞, the coupling is Hilbert–Schmidt, and K(x, x) and K (x, x) are finite at every state, so the standing regularity condition of Theorem H.1 holds. For the forward metric, (16) with (i) and (ii) gives 2

Z

dS (x, y) =

p(z | x) − p(z | y) pX ′ (z)

2 dz ≤

L2 |X | ρ(x, y)2 . c

For the backward metric, the reversed Rkernel at the future state x has density gx (w) = p(x | w)/pX ′ (x) with respect to µ, and |pX ′ (x) − pX ′ (y)| ≤ |p(x | w) − p(y | w)| dµ(w) ≤ Lρ(x, y), so for every w |gx (w) − gy (w)| ≤

1 |pX ′ (y) − pX ′ (x)| C |p(x | w) − p(y | w)| + p(y | w) ≤ L + 2 ρ(x, y), pX ′ (x) pX ′ (x) pX ′ (y) c c

R whence dS ∗ (x, y)2 = (gx − gy )2 dµ ≤ L2 (1/c + C/c2 )2 ρ(x, y)2 . The two bounds make the joint embedding J : (X , ρ) → (J(X ), dJ ) Lipschitz and surjective, and Lipschitz maps do not increase upper box dimension [Falconer, 2014], so d∗ = dimbox (J(X ), dJ ) ≤ dimbox (X , ρ) = d. The hypotheses hold, for example, for a nondegenerate diffusion observed at a fixed lag on a compact manifold: the transition density is then smooth and strictly positive, hence bounded and Lipschitz in each argument.

I

Proof of Theorem 4.1

Every object this proof consumes is constructed in Section H: the forward and backward predictive metrics dS and dS ∗ , the joint embedding J =
(Jf , Jb ) and the joint metric dJ it carries as an isometry, and the intrinsic dimension d∗ = dimbox J(X ), dJ of Theorem H.1. This section takes them as given and proves the budget statement about them. Sketch. Finite box dimension forces total boundedness, so J(X ) is a precompact subset of the Hilbert space ℓ2 ⊕ ℓ2 of box dimension d∗ , and passing to the closure gives a compact set of the same dimension. For any integer m > 2d∗ , embedology [Hunt and Kaloshin, 1999] — the infinite-dimensional form of Sauer et al. [1991] — makes a prevalent, hence at least one, bounded linear π : ℓ2 ⊕ ℓ2 → Rm injective on J(X ) with continuous inverse. This is an existence statement, not an explicit map. The composite ϕ = π ◦ J is dJ -Lipschitz with constant ∥π∥, because J is a dJ -isometry. Injectivity recovers J from ϕ, so every σj ϕj and σj ψj is a function of ϕ: the left singular functions become functions of ϕ(X) and the right ones of ϕ(X ′ ), whence ∥Tϕ ∥HS = ∥T ∥HS . The full argument follows.
Write E = ℓ2 ⊕ℓ2 and recall the joint embedding J(x) = (σj ϕj (x))j , (σj ψj (x))j ∈ E, which is an isometry 2 of (X , dJ ) onto J(X ): ∥J(x) − J(y)∥E = dS (x, y)2 + dS ∗ (x, y)2 = dJ (x, y)2 . The singular system enters only
as a choice of coordinates here. By (23) the same isometry is x 7→ r(x, ·), r∗ (·, x) ∈ L2 (µ′ ) ⊕ L2 (µ), built from the conditional laws alone, and (σj ϕj (x))j is its coordinate vector in the singular basis. Everything below may be read in either presentation.3 By hypothesis d∗ = dimbox (J(X ), dJ ) < ∞, and a set of finite box dimension admits a finite ε-cover for every ε > 0, hence is totally bounded. Since E is complete, J(X ) is 3

The argument can also be closed without the singular system. Since π is injective on the compact
M with continuous inverse and J is an isometry, ϕ is injective on X /{dJ = 0} with Borel inverse, so r(x, x′ ) = H ϕ(x), ϕ(x′ ) for a Borel H. The pair’s 2 density ratio is then σ(ϕ(X)) ⊗ σ(ϕ(X ′ ))-measurable, whence χ2 (ϕ(X); ϕ(X ′ )) = χ2 (X; X ′ ), hence Tϕ HS = ∥T ∥2HS by (13).

28

therefore precompact, and we may replace it by its closure M ⊆ E, a compact set of the same box dimension. (Hilbert–Schmidt decay alone would not give this: the ϕj are orthonormal in L2 (µ), not uniformly bounded, 2 so ∥J(x)∥ = χ2 (p(· | x)∥µ′ ) + χ2 (q(· | x)∥µ) need not be bounded over x.) Fix an integer m > 2d∗ . By the embedding theorem for compact sets of finite box dimension [Sauer et al., 1991, Hunt and Kaloshin, 1999], a prevalent — hence nonempty — set of bounded linear maps π : E → Rm is injective on M with Hölder-continuous inverse on π(M ). Fix such a π and put ϕ = π ◦ J. Admissibility. For x, y ∈ X ,
|ϕ(x) − ϕ(y)| = π J(x) − J(y) ≤ ∥π∥ ∥J(x) − J(y)∥E = ∥π∥ dJ (x, y), so ϕ is dJ -Lipschitz, in particular continuous in the joint metric. ϕ resolves the singular system. Since π is injective on M ⊇ J(X ), J(x) = (π|M )−1 (ϕ(x)), so every coordinate of J — each σj ϕj and each σj ψj — is a continuous function of ϕ. For σj > 0 this makes ϕj measurable with respect to σ(ϕ(X)), that is ϕj ∈ Vϕ , and likewise ψj ∈ Vϕ′ . P The full norm is attained. Writing the singular decomposition T = j σj ⟨ψj , ·⟩L2 (µ′ ) ϕj and using PVϕ ϕj = ϕj , PVϕ′ ψj = ψj for every j with σj > 0, PVϕ T PVϕ′ =

X

D E X σj PVϕ′ ψj , · PVϕ ϕj = σj ⟨ψj , ·⟩ ϕj = T,

j

so ∥Tϕ ∥HS = ∥T ∥HS . Conversely

PVϕ T PVϕ′

j

HS

≤ ∥T ∥HS for every ϕ, since orthogonal projections are 2

2

contractions, so this ϕ attains the largest value the criterion can take, ∥Tϕ ∥HS = ∥T ∥HS , which is (6). In particular ⌊2d∗ ⌋ + 1 coordinates suffice.

J

Experimental Details and Reproducibility

Compute. All experiments run CPU-only on a single Apple M1 laptop (macOS), in Python 3.11 with PyTorch [Paszke et al., 2019]. Library versions are pinned in requirements.txt. Seeds. Data seeds are fixed integer literals: 1 for the ϕ training pairs, 2 for the model-selection split, and 10/11/12 for the readout train/validation/test pools. Model-initialization seeds are derived formulas, 311 + 7i over the ϕ restarts and 17 + 991 r over the VAMPnet restarts, and the warp maps carry their own fixed seeds (7, 11, 99). The dataset is therefore identical across what we report as “seeds”, so the seed spread measures optimizer-initialization variability and not sampling variability. The sole sampling-side error bar in the paper is the resampled label draws of Section 5.5. No torch determinism flags are set. Results are float32 on CPU and machine-stable, but not guaranteed bitwise-portable across BLAS builds. Metrics. Every recovery number is a per-component test R2 : writing y for one target coordinate over the held-out P test pool, Pŷ for its prediction from the learned features, and ȳ for the mean of y on that pool, R2 = 1 − (y − ŷ)2 / (y − ȳ)2 , computed per coordinate and averaged within a component (two coordinates for a cool, one for a bit). It is not clipped, so a representation that predicts a component worse than that component’s own mean scores negative, and a component the representation has dropped reads 0 up to readout noise. Span-side selection uses VAMP-E [Wu and Noé, 2020], which at a rank-k feature pair takes the value 2 tr(U ⊤ C01 V ) − tr(U ⊤ C00 U V ⊤ C11 V ) with the whitening transforms U, V fitted on the training split and the feature covariances — C00 and C11 instantaneous on the present and future features, C01 across the lag — formed on a disjoint one, so it is an unbiased estimate of a population score that overfitting can only depress. Figure S3 reports the fraction built of a singular mode, R2 /σ 2 , where R2 is that of a linear readout of the mode’s future value from the present features and σ 2 is the population R2 available for a unit-variance singular function of singular value σ: one means a fully represented mode, zero an absent one. Hyperparameter searches. Beyond the settings of record, the following ranges were tried during development. Span: trunk widths 642 , 1282 , 2562 , 5122 at ranks k ∈ {12, 30, 50, 80, 100, 120}. Algebra: embedding dimension m ∈ {1, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20}, critic learning rate 2 · 10−4 to 2 · 10−3 in four steps, with 29

held-out Jˆ increasing across the range so the largest was kept, and critic weight decay 0 or 10−4 , the latter never selected. On the warped observation we swept the ϕ learning rate over 10−3 , 3 · 10−4 , 10−4 both cold (held-out Jˆ of 93, 116, 86) and warm-started (140, 161, 173), keeping 3 · 10−4 , and the ϕ trunk over 1282 , 2562 , 5122 , 1283 , 1284 , every widening or deepening of which lowered Jˆ (93, 79, 62, 75, 60), so the warped ϕ stays at 1282 . The readout tunes learning rate {3 · 10−3 , 10−2 } against weight decay {10−5 , 10−3 } on its own validation split. Restarts and stopping epochs were selected on held-out VAMP-E for the span and on held-out Jˆ for the algebra, and no hyperparameter anywhere was selected on a test readout.

30

Related documents

• The Relative Zero Operator: A Causal Nonlocal Operator with Memory in Functional and Spectral Analysis#84639
• Diffusion learning reveals viable parameter manifolds and compensation geometry in biological dynamical systems#1069851
• Input Space Bifurcation Manifolds of Recurrent Neural Networks#492948
• Sparse Representations of Dynamical Networks: A Coprime Factorization Approach#670251
• Sparse Representations of Dynamical Networks: A Coprime Factorization Approach#668186