KPFlow: Operator-Theoretic Gradient Descent

• KPFlow is an operator-theoretic framework that decomposes gradient descent in recurrent systems into dynamic propagation and parameter-filtering components.
• It induces hidden-state updates via the formula δz = -P K P* Err, linking adjoint backpropagation with linearized dynamical propagation.
• The framework enables analysis of multitask alignment and stability in RNNs, and its matrix-free PyTorch implementation handles large-scale trajectory data.

KPFlow denotes an operator-theoretic framework for gradient-descent learning in recurrent dynamical systems, the pair of operators $K$ and $P$ that define that framework, and a PyTorch-based software package implementing the resulting diagnostics in matrix-free form (Hazelden et al., 8 Jul 2025). In the formulation introduced in “KPFlow: An Operator Perspective on Dynamic Collapse Under Gradient Descent Training of Recurrent Networks” (Hazelden et al., 8 Jul 2025), a gradient step induces a hidden-state update

$\delta z = - P K P^*\, Err,$

so that learning is decomposed into adjoint backpropagation through the linearized dynamics, parameter-space filtering, and forward propagation through state space. The framework is intended for finite, nonlinear recurrent models, including RNNs, GRUs, Neural ODEs, and, in the appendix, more general dynamical systems such as Hodgkin–Huxley networks.

1. Definition and problem setting

KPFlow is introduced simultaneously as a theoretical decomposition, an operator pair, and a software package (Hazelden et al., 8 Jul 2025). Its central domain is recurrent dynamical systems trained by gradient descent, written in continuous time as

$\frac{d}{dt} z(t|x) = f(z(t|x), x(t)\mid \theta),$

with supervised trajectory loss

$$L(\theta) := \mathbb{E}_{x \sim X}\Big[\int_0^{t_{end}} \ell(z(t|x,\theta), y^*(t|x))\, dt\Big].$$

The framework is motivated by three recurrent-learning phenomena emphasized in the paper: representation dynamics, latent collapse or simplicity bias, and multitask alignment or interference.

The stated claim is that these phenomena cannot be understood from task structure alone. The architecture and parameterization induce a learning operator that constrains which hidden-state perturbations gradient descent can realize. In this setting, KPFlow separates two influences: the geometry of the parameterization and the propagation structure of the recurrent dynamics. This separation is the basis for both the theoretical analysis and the accompanying diagnostics.

The formal state space is a trajectory space $T$. For a fixed trial $x$,

$T_x = L^2([0,t_{end}]; \mathbb{R}^n),$

and over all trials $T$ is the direct sum or direct integral of these spaces. After discretization, trajectories are represented as tensors of shape $[B,T,H]$, with batch, time, and hidden dimensions.

2. Gradient-flow factorization

The derivation begins from gradient-based training via adjoint equations (Hazelden et al., 8 Jul 2025). For a trial $P$0, the paper defines the adjoint

$P$1

and the instantaneous loss sensitivity

$P$2

For MSE with linear readout, the appendix gives

$P$3

The adjoint evolves backward according to

$P$4

A running parameter gradient $P$5 is defined, with

$P$6

To propagate perturbations in state space, the paper introduces the state-transition matrix

$P$7

where $P$8. This matrix supplies the linearized recurrent propagator that connects KPFlow to Lyapunov stability analysis.

From these ingredients, the paper derives the hidden-state perturbation induced by one gradient step and factors it into a backward propagator, a parameter kernel, and a forward propagator. Proposition 1 states the resulting decomposition: $P$9 This is the defining identity of KPFlow.

The paper also situates this factorization among special cases. For feedforward or no time-stepping settings, $\delta z = - P K P^*\, Err,$0; for optimize-then-discretize or direct perturbation of flow, $\delta z = - P K P^*\, Err,$1; and for standard deep-learning discretize-then-optimize, $\delta z = - P K P^*\, Err,$2. This places KPFlow as a recurrent generalization of an NTK-style perspective.

3. The operators $\delta z = - P K P^*\, Err,$3 and $\delta z = - P K P^*\, Err,$4

The propagation operator $\delta z = - P K P^*\, Err,$5 is defined by

$\delta z = - P K P^*\, Err,$6

Its implementation-oriented form is the variational equation

$\delta z = - P K P^*\, Err,$7

with $\delta z = - P K P^*\, Err,$8. Operationally, if $\delta z = - P K P^*\, Err,$9 is a perturbation to the hidden-state forcing term, $\frac{d}{dt} z(t|x) = f(z(t|x), x(t)\mid \theta),$0 is the resulting perturbation in the hidden trajectory.

The parameter operator $\frac{d}{dt} z(t|x) = f(z(t|x), x(t)\mid \theta),$1 is

$\frac{d}{dt} z(t|x) = f(z(t|x), x(t)\mid \theta),$2

where $\frac{d}{dt} z(t|x) = f(z(t|x), x(t)\mid \theta),$3. The integral-plus-expectation first produces a parameter-space perturbation; $\frac{d}{dt} z(t|x) = f(z(t|x), x(t)\mid \theta),$4 then maps that perturbation back into hidden-state forcing space. In the language of the paper, $\frac{d}{dt} z(t|x) = f(z(t|x), x(t)\mid \theta),$5 is the operator induced by the parameterization of the dynamics.

The paper identifies $\frac{d}{dt} z(t|x) = f(z(t|x), x(t)\mid \theta),$6 as the recurrent or dynamical analogue of a Neural Tangent Kernel-like operator, but on trajectory space rather than input-output space (Hazelden et al., 8 Jul 2025). In discretized form,

$\frac{d}{dt} z(t|x) = f(z(t|x), x(t)\mid \theta),$7

Hence each kernel block compares parameter sensitivities across trials, times, and hidden coordinates, and the singular vectors or eigenfunctions of $\frac{d}{dt} z(t|x) = f(z(t|x), x(t)\mid \theta),$8 are themselves trajectory-shaped objects.

Theorem 1 further decomposes the parameter operator across parameters or computation-graph views: $\frac{d}{dt} z(t|x) = f(z(t|x), x(t)\mid \theta),$9 Each $$L(\theta) := \mathbb{E}_{x \sim X}\Big[\int_0^{t_{end}} \ell(z(t|x,\theta), y^*(t|x))\, dt\Big].$$0 is a kernel-induced Hilbert–Schmidt operator and is therefore compact, continuous, self-adjoint, and positive definite. Its effective rank is tied to the effective dimension of the dynamical signal $$L(\theta) := \mathbb{E}_{x \sim X}\Big[\int_0^{t_{end}} \ell(z(t|x,\theta), y^*(t|x))\, dt\Big].$$1 to which weight $$L(\theta) := \mathbb{E}_{x \sim X}\Big[\int_0^{t_{end}} \ell(z(t|x,\theta), y^*(t|x))\, dt\Big].$$2 is applied. This is the basis of the paper’s claim that low-dimensionality can arise from network structure before task structure is considered.

Theorem 2 gives an analogous structural result for the propagator: $$L(\theta) := \mathbb{E}_{x \sim X}\Big[\int_0^{t_{end}} \ell(z(t|x,\theta), y^*(t|x))\, dt\Big].$$3 where $$L(\theta) := \mathbb{E}_{x \sim X}\Big[\int_0^{t_{end}} \ell(z(t|x,\theta), y^*(t|x))\, dt\Big].$$4 is induced by the fundamental matrix $$L(\theta) := \mathbb{E}_{x \sim X}\Big[\int_0^{t_{end}} \ell(z(t|x,\theta), y^*(t|x))\, dt\Big].$$5 and $$L(\theta) := \mathbb{E}_{x \sim X}\Big[\int_0^{t_{end}} \ell(z(t|x,\theta), y^*(t|x))\, dt\Big].$$6 is a Volterra integral operator. The proof writes

$$L(\theta) := \mathbb{E}_{x \sim X}\Big[\int_0^{t_{end}} \ell(z(t|x,\theta), y^*(t|x))\, dt\Big].$$7

This ties $$L(\theta) := \mathbb{E}_{x \sim X}\Big[\int_0^{t_{end}} \ell(z(t|x,\theta), y^*(t|x))\, dt\Big].$$8 directly to the linearized dynamics, Lyapunov directions, and finite-time growth or decay structure.

4. Mechanism of representation formation and collapse

The mechanistic interpretation proposed in the paper proceeds in four stages (Hazelden et al., 8 Jul 2025). First, task error generates an error trajectory $$L(\theta) := \mathbb{E}_{x \sim X}\Big[\int_0^{t_{end}} \ell(z(t|x,\theta), y^*(t|x))\, dt\Big].$$9. Second, backward propagation through $T$0 yields the adjoint signal $T$1, identifying hidden-state directions whose perturbation would reduce loss. Third, filtering through $T$2 restricts those directions to the subset realizable through parameter changes. Fourth, forward propagation through $T$3 converts those realizable perturbations into hidden-state updates $T$4.

In this account, learned representation is governed not merely by the target task but by the geometry of

$T$5

The paper argues that latent collapse and low-dimensional dynamical motifs arise because $T$6 amplifies or suppresses perturbations according to recurrent stability structure, while $T$7 is often effectively low-rank because parameter changes induce only a restricted family of trajectory perturbations. This suggests that gradient descent preferentially follows dominant singular or eigen-directions of the composite operator, producing low-dimensional attractors, shared subspaces across conditions, and a recurrent simplicity bias.

The appendix contrasts parameterized learning with direct perturbation of the flow. If one could perturb the flow directly, the correction would be

$T$8

Under parameterization, however,

$T$9

so the effective correction directions need not coincide with those of $x$0. The paper calls this distortion gradient misdirection: parameterization rotates or distorts the ideal correction dictated by the adjoint and the dynamics.

Corollary 1 extends the same structure to outputs. For MSE with linear readout, the output update is

$x$1

Thus the operator decomposition governs both latent-state motion and output correction.

5. Multitask alignment and interference

A second major application of KPFlow is the analysis of multitask recurrent training (Hazelden et al., 8 Jul 2025). The paper proposes task- or subtask-specific operators such as $x$2, intended to quantify how the parameter structure associated with one task interacts with the hidden-state geometry of another.

The formal alignment diagnostic is based on the Rayleigh quotient

$x$3

with inner product over time, trials, and hidden units. Large $x$4 indicates strong alignment of $x$5 with dominant eigendirections of $x$6. For multitask analysis, the appendix proposes matrices of the form

$x$7

where $x$8 is the subtask error and $x$9 its adjoint signal. This yields an interference or alignment matrix that quantifies constructive sharing versus conflict.

The paper also reports empirical alignment measurements on the resulting hidden-state objectives $T_x = L^2([0,t_{end}]; \mathbb{R}^n),$0, using cumulative cosine similarities across subtask pairs. The stated point is that KPFlow can explain why some subtasks share attractors or common subspaces, and why others interfere constructively or destructively depending on initialization scale and recurrent dynamics.

In GRU multitask experiments, small initial recurrent scale $T_x = L^2([0,t_{end}]; \mathbb{R}^n),$1 leads to stronger alignment into particular organizations, notably Pro/Anti sharing, whereas large $T_x = L^2([0,t_{end}]; \mathbb{R}^n),$2 yields less aligned and more distributed updates early in training. These findings are presented as empirical support for the operator analysis rather than as formal theorems.

6. Software implementation and empirical validation

The KPFlow package implements the operator framework in PyTorch and is designed around matrix-free operator action rather than explicit matrix construction (Hazelden et al., 8 Jul 2025). For discretized trajectories of shape $T_x = L^2([0,t_{end}]; \mathbb{R}^n),$3, explicit operator matrices would have size $T_x = L^2([0,t_{end}]; \mathbb{R}^n),$4, so the package applies operators implicitly.

For $T_x = L^2([0,t_{end}]; \mathbb{R}^n),$5, the appendix describes a two-stage procedure. First compute

$T_x = L^2([0,t_{end}]; \mathbb{R}^n),$6

then evaluate

$T_x = L^2([0,t_{end}]; \mathbb{R}^n),$7

This is implemented with vjp for the first step and jvp for the second. For $T_x = L^2([0,t_{end}]; \mathbb{R}^n),$8, the package provides two routes: a linear approach that solves the variational equation, and a direct approach based on

$T_x = L^2([0,t_{end}]; \mathbb{R}^n),$9

valid for sufficiently small perturbations. A QR-stabilized method based on Lyapunov-exponent computation is also implemented for more stable propagation through Jacobian products.

The operators are wrapped in a generic Operator class that can be converted to a SciPy LinearOperator, enabling matrix-free SVD, eigendecomposition, and related spectral analysis. The package also includes an AveragedOperator class for averaging over selected axes, such as trials, times, or hidden units, in order to compute consensus spectra.

The experimental program in the appendix covers tanh RNNs, GRUs, RNNs with ReLU, Neural ODE-compatible analysis, and a derivation for Hodgkin–Huxley networks. Tasks include Memory-Pro and multitask settings combining Memory-Pro, Memory-Anti, Delay-Pro, and Delay-Anti. A typical setup uses hidden size $T$0, batch trials, and ADAM in experiments, with metrics including convergence time, effective rank of $T$1, hidden-state dimensionality, Lyapunov dimension, cosine similarity, and operator alignment.

The reported findings are consistent across several regimes. First, the effective rank of $T$2 is often much lower than that of $T$3, so parameterization acts as a bottleneck. Second, larger initial recurrent scale $T$4 tends to increase dynamical richness and can speed learning up to a point, with low $T$5 associated with more regular and lower-dimensional dynamics and high $T$6 with more chaotic and higher-dimensional dynamics. Third, the Rayleigh quotient

$T$7

typically increases as loss decreases, suggesting growing alignment between error signals and learnable directions of $T$8. Fourth, ReLU RNNs exhibit less bottlenecking than tanh RNNs, with more chaotic and higher-dimensional dynamics and larger effective rank of $T$9. These are presented as empirical findings rather than proofs.

7. Scope, limitations, and related nomenclature

The framework is derived for gradient flow or infinitesimal GD-style updates, whereas the experiments use ADAM, which the paper notes is related but not identical (Hazelden et al., 8 Jul 2025). The operator $[B,T,H]$0 is fundamentally linearized, and the direct perturbation approximation

$[B,T,H]$1

is valid only for sufficiently small perturbations. All numerical analysis is performed on discretized trajectories even though the conceptual framework is continuous-time.

The paper also notes that Theorem 1 gives exact equality between effective rank and dynamical-signal dimension only in certain cases, while more general settings yield upper bounds. Some multitask conclusions remain empirical and may depend on initialization regime. The exact relation between the SVD of $[B,T,H]$2 and classical Lyapunov spectra is identified as future work.

A later paper, “The Global Empirical NTK: Self-Referential Bias and Dimensionality of Gradient Descent Learning” (Hazelden et al., 9 May 2026), introduced kpflow as a matrix-free library for analyzing the Global Empirical NTK with the related factorization

$[B,T,H]$3

That work extends the operator vocabulary to a global-state NTK setting including RNNs and transformers. This suggests continuity in notation and operator structure, but the 2025 recurrent-learning paper defines KPFlow specifically through the hidden-state update

$[B,T,H]$4

and the corresponding PyTorch analysis package.

The term should also be distinguished from unrelated or differently scoped names in the literature, including Kernel Flows-based kernel tuning for MSPC or K-PLS (Duma et al., 2 May 2025, Duma et al., 2023), the generative method “K-Flow” (Du et al., 27 Apr 2025), and the C++ task-parallel framework “Pipeflow” (Chiu et al., 2022). Within its own primary usage, KPFlow denotes an operator perspective on how gradient descent shapes recurrent representations, with the central claim that latent collapse, task sharing, and interference are properties of the learned task together with the architecture-induced operator geometry, not of the task alone (Hazelden et al., 8 Jul 2025).