Paper deep dive
MaRK: Markov-adapted Recurrent Kernels for Dynamic Operator Conditioning in State Space Models
Syed Ibrahim Omer, Ginny Y. Wong, Xiangyu Zhao
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 10/9/2026, 7:32:28 AM
Summary
The paper introduces MaRK (Markov-adapted Recurrent Kernels), a framework for dynamic operator conditioning in State Space Models (SSMs). MaRK modulates the internal recurrence parameters (A, B, C, D, Δ) of a frozen SSM backbone (Hydra) using context vectors derived from diffusion timesteps. It employs three adapter geometries (Hypernet, Chebyshev, DCT) to map context to bounded parameter shifts, ensuring stability via Affine Quadratic Stability (AQS) certificates. The method enables parameter-efficient fine-tuning for iterative diffusion generation, with the Chebyshev variant achieving the best performance.
Entities (9)
Relation Signals (9)
MaRK → modulates → State Space Models
confidence 95% · MaRK (Markov-adapted Recurrent Kernels), a dynamic operator-conditioning framework that maps context vectors directly into bounded modulations of a frozen SSM's recurrence
MaRK → instantiatedon → Hydra
confidence 92% · We instantiate MaRK on a frozen 111M-parameter Hydra SSM backbone
MaRK → usesadapter → Hypernet
confidence 90% · study three adapter geometries: Hypernet, Chebyshev polynomial, and Discrete Cosine Transform kernels.
MaRK → usesadapter → DCT
confidence 90% · study three adapter geometries: Hypernet, Chebyshev polynomial, and Discrete Cosine Transform kernels.
MaRK → usesadapter → Chebyshev
confidence 90% · study three adapter geometries: Hypernet, Chebyshev polynomial, and Discrete Cosine Transform kernels.
Chebyshev → achievesbestperformance → MaRK
confidence 88% · Empirically, the Chebyshev variant yields the strongest performance, achieving an average validation loss of 2.55
MaRK → providesstability → Affine Quadratic Stability
confidence 88% · The bounded recurrence parameterization further yields an analytic Affine Quadratic Stability certificate for the modulated recurrence.
MaRK → basedontheory →
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:State Space Models (SSMs) offer an efficient alternative to Transformers for sequence modeling, yet conditioning pre-trained SSMs for iterative generation typically operates outside the recurrent operator, through input injection or activation modulation. While such mechanisms expose the model to conditioning information, they leave the underlying temporal dynamics fixed. We introduce MaRK (Markov-adapted Recurrent Kernels), a dynamic operator-conditioning framework that maps context vectors directly into bounded modulations of a frozen SSM's recurrence ($A$), read-in ($B$), read-out ($C$), skip ($D$), and discretization ($\Delta$) parameters. Viewed through the lens of LPV-SSM systems, MaRK induces a context-indexed family of Markov parameter sequences, allowing each diffusion timestep to reshape the model's input-output memory kernel. We instantiate MaRK on a frozen 111M-parameter Hydra SSM backbone and study three adapter geometries: Hypernet, Chebyshev polynomial, and Discrete Cosine Transform kernels. Since these adapters modify the Markov parameter sequence through low-rank auxiliary maps on the frozen backbone, parameter-efficient fine-tuning arises as a structural consequence of the adaptation mechanism itself, requiring only 6.3--11M trainable auxiliary parameters to transition from a bidirectional objective to an iterative diffusion regime. The bounded recurrence parameterization further yields an analytic Affine Quadratic Stability certificate for the modulated recurrence. Through synthetic LPV recovery experiments and Markov-operator diagnostics, we show that MaRK recovers coordinate-invariant temporal operators under matched assumptions and produces distinct, stable timestep-conditioned memory profiles. Empirically, the Chebyshev variant yields the strongest performance, achieving an average validation loss of 2.55, followed by the DCT (2.59) and Hypernet (3.77) geometries.
Tags
Links
- Source: https://arxiv.org/abs/2610.09092v2
- Canonical: https://arxiv.org/abs/2610.09092v2
Trouble viewing inline? Open PDF directly →
Full Text
106,655 characters extracted from source content.
Expand or collapse full text
MaRK: Markov-adapted Recurrent Kernels for Dynamic Operator Conditioning in State Space Models Syed Ibrahim Omer †thanks: Corresponding Author: ibrahsyed2-c@my.cityu.edu.hk Affiliation: Department of Data Science Affiliation: City University of Hong Kong Ginny Y. Wong Affiliation: NVIDIA AI Technology Center Affiliation: NVIDIA, Santa Clara, USA Xiangyu Zhao Affiliation: Department of Data Science Affiliation: City University of Hong Kong Abstract State Space Models (SSMs) offer an efficient alternative to Transformers for sequence modeling, yet conditioning pre-trained SSMs for iterative generation typically operates outside the recurrent operator, through input injection or activation modulation. While such mechanisms expose the model to conditioning information, they leave the underlying temporal dynamics fixed. We introduce MaRK (Markov-adapted Recurrent Kernels), a dynamic operator-conditioning framework that maps context vectors directly into bounded modulations of a frozen SSM’s recurrence (A), read-in (B), read-out (C), skip (D), and discretization (Δ ) parameters. Viewed through the lens of Linear Parameter-Varying systems, MaRK induces a context-indexed family of Markov parameter sequences, allowing each diffusion timestep to reshape the model’s input-output memory kernel. We instantiate MaRK on a frozen 111M-parameter Hydra SSM backbone and study three adapter geometries: Hypernet, Chebyshev polynomial, and Discrete Cosine Transform kernels. Since these adapters modify the Markov parameter sequence through low-rank auxiliary maps on the frozen backbone, parameter-efficient fine-tuning arises as a structural consequence of the adaptation mechanism itself, requiring only 6.3–11M trainable auxiliary parameters to transition from a bidirectional objective to an iterative diffusion regime. The bounded recurrence parameterization further yields an analytic Affine Quadratic Stability certificate for the modulated recurrence. Through synthetic LPV recovery experiments and Markov-operator diagnostics, we show that MaRK recovers coordinate-invariant temporal operators under matched assumptions and produces distinct, stable timestep-conditioned memory profiles. Empirically, the Chebyshev variant yields the strongest performance, achieving an average validation loss of 2.55, followed by the DCT (2.59) and Hypernet (3.77) geometries. Together, these results provide initial evidence that dynamic operator modulation is a principled operator-level conditioning mechanism for adapting SSMs beyond input-stream injection and adaptive normalization. 1 Introduction State Space Models (SSMs) have emerged as highly efficient, continuous-time alternatives to Transformers for large-scale sequence modeling [1, 2, 3] . However, adapting sequence models for conditional generation tasks (e.g., discrete diffusion) traditionally relies on post-hoc interventions [4, 5, 6, 7, 8]. This challenge is commonly addressed by two external conditioning paradigms: Input-Stream Injection [4] and Adaptive Layer Normalization (AdaLN) [9, 10, 11]. Input-stream injection concatenates conditioning vectors directly into the sequence; while computationally cheap, it demands that the state-space dynamics route this context across the entire sequence length. For continuous-time SSMs, this relies on learned hidden states carrying context forward, which can dilute the prompt signal and induce temporal lag. Conversely, AdaLN dynamically modulates structural layer activations. While effective for Transformers, applied to SSMs, AdaLN operates externally to the core continuous-time sequence operator. It scales layer inputs and outputs but does not explicitly parameterize the underlying temporal dynamics (the A, B, C, and Δ matrices) as functions of the timestep context. Both paradigms condition the computation without directly parameterizing the Markov parameter sequence (ℋ(ct)H(c_t)), the complete characterization of an SSM’s input-output temporal behavior. To bridge this gap, we introduce MaRK (Markov-adapted Recurrent Kernels), a structurally principled operator-level conditioning approach: Dynamic Operator Modulation. Rather than relying only on external activation normalization or input additions, MaRK moves the conditioning map into the temporal recurrence parameters themselves. By projecting the context vector into bounded multi-dimensional shifts of the fundamental continuous-time matrices, it dynamically modulates the internal recurrence (A), read-in (B), read-out (C), skip connections (D), and discretization steps (Δ ) of a frozen base model [12]. We apply our method to the Hydra BERT SSM model, a bi-directional extension of the Mamba model parameterized as a quasi-separable matrix mixer [13]. Through MaRK, the continuous-time matrices of the pre-trained Hydra BERT SSM depend on external conditioning signals (diffusion timesteps), allowing the system’s temporal properties to vary dynamically at each step. Because MaRK relies on bounded continuous modulations via low-rank geometries, the same design yields a parameter-efficient adaptation route in our Hydra instantiation, providing a dynamic alternative to static low-rank adaptation methods [14]. We present three kernel geometries: Hypernet, Chebyshev Polynomial, and Discrete Cosine Transform (DCT) [15] kernels, that map conditioning contexts to scalar shifts and low-rank scales. By bounding parameter modulations via hyperbolic tangent saturation and softplus constraints, MaRK admits Affine Quadratic Stability (AQS) certificates for the modulated recurrence [16]. Applied to a 2323-layer Hydra SSM, our adapters use 6.3–11M trainable parameters to condition the frozen 111M base structure, turning a computationally expensive pre-training problem into a targeted fine-tuning regime [13]. We evaluate this direct conditioning approach within an iterative mask-based diffusion framework using stability certificates, synthetic operator recovery, Markov-operator diagnostics, and internal adapter ablations. MaRK specifically adapts the frozen Hydra backbone to an iterative diffusion process, employing a Context-Adaptive Token-Level Noise Rescheduling (CART) objective to provide token-level, context-aware training signal within our iterative refinement framework [8]. Contributions. (i) LPV-SSM Operator Conditioning: We cast conditional generation as LPV control over SSM operators, directly modulating A, B, C, D, and Δ as an operator-level perspective beyond input injection and AdaLN. (i) MaRK Kernel Family: We design three structural adapter geometries (Hypernet, Chebyshev, and DCT) that map conditioning vectors to low-rank modulations. (i) Stability Certificates and Operator Recovery: We use bounded parameterization to obtain analytic AQS certificates and study Markov-operator recovery under matched synthetic assumptions. (iv) Markov Operator Adaptation: We provide evidence that MaRK conditions the Markov parameter sequence, with timestep-specific lag-norm profiles revealing distinct temporal-memory signatures. 2 Linear Parameter Varying (LPV) Foundation We adapt the base time-invariant State Space Model using the theory of Linear Parameter Varying (LPV) models [16, 17]. Therefore, we view each modulated SSM layer as a continuous-time state-space system: x˙(τ) x(τ) =A(t)x(τ)+B(t)u(τ) =A(c_t)\,x(τ)+B(c_t)\,u(τ) (1) y(τ) y(τ) =C(t)x(τ)+D(t)u(τ) =C(c_t)\,x(τ)+D(c_t)\,u(τ) (2) where x(τ)∈ℝnx(τ) ^n is the state, u(τ)∈ℝmu(τ) ^m the input, y(τ)∈ℝmy(τ) ^m the output, n is the state dimension, m the embedding dimension, tc_t is a dynamic conditioning vector tied to the diffusion timestep t obtained using standard DDPM/DiT sinusoidal (Fourier) encoding and two-layer MLP [4, 9], x˙(τ) x(τ) is the derivative of the state x(τ)x(τ) with respect to continuous internal-time τ (i.e. dxdτ dxdτ). MaRK modulates the continuous-time parameters Q∈A,B,C,D,ΔQ∈\A,B,C,D, \ using learned mapping functions ψscale(⋅) _scale(·) and ψshift(⋅) _shift(·). We use the primed notation to denote the modulated parameters Q′∈A′,B′,C′,D′,Δ′Q ∈\A ,B ,C ,D , \. For diagonalized scalar parameters Θ∈Alog,Δbias,D ∈\A_log, _bias,D\ we apply bounded shifts to maintain stability, where βΘ _ denotes the scalar magnitude of the bounded shifts: Θ′ =Θ+βΘtanh(ψshift(t)) = + _ ( _shift(c_t) ) (3) For dense projection matrices Φ∈B,C ∈\B,C\ we apply bounded scales and shifts, where αΦ _ , βΦ _ denote the scalar magnitude of the bounded scales and shifts respectively. To avoid the curse of dimensionality for large matrices B,C∈ℝdB,C ^d (where d=ngroups×dstated=n_groups× d_state; ngroupsn_groups is the number of SSM groups and dstated_state is the state dimension), modulation is applied via a rank-r factorized projection to these matrices: Φ′ =Φ⊙[1+αΦtanh(ψscale(t))]+βΦtanh(ψshift(t)) = [1+ _ ( _scale(c_t) ) ]+ _ ( _shift(c_t) ) (4) Remark (Operator vs. Activation Control). While §1 contrasts these paradigms conceptually, the structural distinction is formalized in Eq. (3)–(4). By modulating the infinitesimal generator A and discretization Δ rather than applying scalar affine transforms (AdaLN), MaRK directly parameterizes the system’s eigenvalue distribution and frequency response. This makes the context tc_t explicit in the temporal kernel’s decay properties, rather than treating it only as an external signal bias. Finally, we obtain the LPV-SSM form from equations (1) and (2): x˙(τ) x(τ) =A′¯(t)x(τ)+B′¯(t)u(τ) = A (c_t)\,x(τ)+ B (c_t)\,u(τ) (5) y(τ) y(τ) =C′(t)x(τ)+D′(t)u(τ) =C (c_t)\,x(τ)+D (c_t)\,u(τ) (6) where A′¯ A and B′¯ B are the discretized versions of the continuous-time parameters A′A and B′B , obtained via the Zero-Order Hold (ZOH) [12]. The complete input-output behavior of this LPV-SSM is fully characterized by the Markov parameter, where ℋ(t)H(c_t) is the Markov Parameter Sequence, hk(τ)h_k(τ) is the kthk^th Markov Parameter, and k is the shift index [18]: ℋ(t)=hk(τ)=C′(t)A′k¯(t)B′¯(t)k=0∞ (c_t)= \h_k(τ)=C (c_t)\, A \,k(c_t)\, B (c_t) \_k=0^∞ (7) This object is the natural invariant of the system: two SSMs with different internal coordinates but identical Markov sequences are behaviorally indistinguishable. Crucially, by making each matrix in this product condition-dependent, MaRK adapts the full temporal operator at each generation step rather than acting externally on activations or inputs. We term this approach Markov-adapted Recurrent Kernels (MaRK). 2.1 Constructive Stability via Log-Space Parameterization. The bounded modulations above support a constructive stability argument through the ZOH (Zero-Order Hold) discretization chain. In Hydra’s diagonal SSM, each recurrence eigenvalue A¯i A_i is produced by the composition: A A =−exp(Alog′)<0 =- (A_ )<0 (8) Δ =softplus(Δinput+Δbias′)>0 =softplus( _input+ _bias )>0 (9) A¯i A_i =exp(Ai⋅Δi)∈(0,1) = (A_i· _i)∈(0,1) (10) Equation (8) follows from exp(x)>0 (x)>0 for all x∈ℝx , hence −exp(Alog′)<0- (A_ )<0. Equation (9) holds because softplus(x)=log(1+ex)>0softplus(x)= (1+e^x)>0 for all x∈ℝx . Finally, Equation (10) follows since Ai<0A_i<0 and Δi>0 _i>0 imply Ai⋅Δi<0A_i· _i<0, and therefore A¯i=exp(Ai⋅Δi)∈(0,1) A_i= (A_i· _i)∈(0,1). Proposition 2.1 (Constructive AQS stability). For the MaRK-modulated diagonal SSM, the identity matrix P=IP=I is a valid common Lyapunov function satisfying the AQS condition A¯(t)⊤PA¯(t)−P⪯−εI A(c_t) P\, A(c_t)-P - I for all admissible conditioning contexts tc_t, with stability margin ε=1−maxiA¯i2>0 =1- _i A_i^2>0. Proof sketch. Since A is diagonal, the Lyapunov condition reduces element-wise to A¯i2−1<0 A_i^2-1<0 for each head i. By Equations (8)–(10), every element A¯i∈(0,1) A_i∈(0,1) for admissible conditioning contexts. Because A¯i(δ) A_i(δ) is monotonically decreasing in the shift δ (larger shift → more negative A → more decay), the worst-case spectral radius occurs at δ=−βAδ=- _A, bounded away from 11 by the finite hyper-rectangle of modulations. ∎ 3 Method 3.1 Pretrained Hydra BERT SSM backbone We instantiate MaRK on top of a pre-trained Hydra BERT backbone, a bidirectional extension of Mamba implemented through a quasiseparable matrix mixer [13]. Hydra adopts the “matrix mixer” view of sequence mixing, in which a sequence mixer is a structured linear map Y=MXY=MX on the token sequence. In contrast to ad-hoc bidirectional SSM heuristics that combine separate forward/backward models by addition or concatenation [19, 20], Hydra parameterizes a single mixer matrix whose strictly lower and upper triangles correspond to forward and backward semiseparable (SSM) structure, while the diagonal remains free and acts as a learned residual pathway. Concretely, Hydra’s quasiseparable mixer admits the form (cf. Eq. (3) in [13]) mij=→i⊤→i:j×→ji>jδi=j←j⊤←j:i×←i<jm_ij= cases c_i \, A_i:j^×\, b_j&i>j\\ _i&i=j\\ c_j \, A_j:i^×\, b_i&i<j cases (11) where δi _i is a scalar, bi,ci∈ℝN×1b_i,c_i ^N× 1, and Ai∈ℝN×NA_i ^N× N. Since this matrix features non-zero upper triangle components, it enables bidirectional processing [13]. In our experiments, we use the released Hydra-BERT checkpoint at BERT scale (111M parameters, 23 layers) [13]. Hydra is a natural backbone for our setting because discrete diffusion language modeling is inherently bidirectional (tokens are refined using both left and right context), yet our objective requires timestep-dependent control over the internal temporal operator. MaRK targets exactly the parameters that define Hydra’s SSM dynamics (the continuous-time recurrence and discretization operators: A,B,C,D, and ΔA,B,C,D, and underlying the quasi-separable factors) and makes them functions of the diffusion context tc_t. By doing so, we demonstrate the ability of MaRK to also function as a fine-tuning method to convert an existing bidirectional objective into a context-aware diffusion objective. Throughout, we treat the Hydra backbone as frozen and learn only the MaRK adapters, isolating the effect of dynamic operator modulation from changes in the base representation. 3.2 Discrete Diffusion Language Modeling We evaluate MaRK in the masked discrete diffusion framework used by recent diffusion LLMs [8, 7, 5]. Let 0∈Lx_0 ^L be a token sequence of length L from the data distribution, and let ℳM denote an absorbing noise token. The forward process independently corrupts tokens according to a timestep t by masking each position with probability t (equivalently, keeping it unmasked with probability α(t)=1−tα(t)=1-t) [7]: q(t∣0)=∏i=1Lq(xt,i∣x0,i),q(xt,i=ℳ∣x0,i)=t,q(xt,i=x0,i∣x0,i)=1−t.q(x_t _0)= _i=1^Lq(x_t,i x_0,i), q(x_t,i=M x_0,i)=t,\;\;q(x_t,i=x_0,i x_0,i)=1-t. (12) The reverse model is a mask predictor, pθ(0∣t,t)p_θ(x_0 _t,t), that predicts all masked positions in parallel using full bidirectional context. Concretely, we parameterize pθp_θ with the Hydra backbone, and provide the diffusion context (timestep embedding) as the conditioning vector tc_t that drives MaRK operator modulation that conditions pθp_θ based on the diffusion context. Following [8, 7], we train by sampling a clean sequence 0x_0, sampling t∼(0,1)t (0,1), then drawing t∼q(t∣0)x_t q(x_t _0), and minimizing the negative log-likelihood on only the masked tokens: ℒdiff(θ)=0,t,t[w(t)∑i=1L[xt,i=ℳ](−logpθ(x0,i∣t,t))],L_diff(θ)=E_x_0,\,t,\,x_t [w(t) _i=1^L1[x_t,i=M]\; (- p_θ(x_0,i _t,t) ) ], (13) where w(t)w(t) is a timestep-dependent reweighting term derived from the diffusion ELBO for particular schedules [7], and can be further adapted at the token level based on local context informativeness [8]. In our implementation, this context-adaptive reweighting is instantiated by CART (§3.3), which emphasizes masked tokens that are surrounded by more usable (unmasked) context. 3.3 Training Objective To optimize training, we utilize Context-Adaptive Token-level Noise Rescheduling (CART), following the token-level rescheduling paradigm introduced in Dream [8]. Rather than using a single global timestep-dependent weight w(t)w(t) for all masked tokens at diffusion time t (Eq. (13)), CART applies a per-token weight wi(t,t)w_i(t,x_t) that depends on how much usable (unmasked) context surrounds token i in the partially-masked sequence tx_t. This matters because two masked tokens at the same t can have very different predictability: a token with many nearby visible neighbors is effectively less ambiguous and should contribute a stronger learning signal, while a token whose neighborhood is mostly masked is inherently underdetermined and over-penalizing it injects gradient noise and encourages collapse toward “easy” tokens. CART operationalizes “how much context” as a distance-weighted sum of the visible tokens: each masked position receives one symmetric geometric kernel centered on every unmasked position, so nearer visible tokens contribute more and the weight grows with the amount of nearby visible context [8]. The resulting token-wise weights multiply the masked-token cross-entropy terms before averaging. We provide the full mathematical definition and its relationship to Eq. (13) in Appendix A.3. 3.4 MaRK Kernel Adapters MaRK proposes the mapping ψ(⋅)ψ(·) that serves as the system’s scheduling function, translating high-dimensional semantic context into the low-dimensional manifold of stable dynamical parameters. By varying the basis of this mapping, we can control the spectral regularity of the adaptation process. The certified setting uses the bounded modulation structure from equations (3) and (4), and the three adapter geometries are ablated head-to-head in §4.4. Hypernet Adapter: The Hypernet maps the hidden state h (obtained from the tc_t via a mapping) directly to parameter shifts using linear projections W, to get the left and right low-rank matrices U and V, functioning as the most parameter-efficient baseline (≈ 6.3M parameters): Θshift _shift =WΘh =W_ h (14) UVΦscale/shift UV_ _scale/shift =WΦscale/shifth =W_ _scale/shifth (15) The bounded scales control the gain amplification avoiding explosive states, whereas the shifts are kept optionally left unconstrained for the Hypernet kernel only, to be treated as a residual correction and act as the most expressive baseline. Chebyshev Polynomial Adapter: To enforce smoother gradients, this kernel uses orthogonal polynomial basis functions evaluated at a latent coordinate z=Wzhz=W_zh. The Chebyshev polynomials Tk(z)T_k(z) (up to degree K) determine the continuous shifts (≈ 11M parameters): Θshift=∑k=1K[WΘ(h)]k⋅Tk(z) _shift= _k=1^K[W_ (h)]_k· T_k(z) (16) For matrix elements, we utilize the basis mean T¯k=1nheads∑iTk(zi) T_k= 1n_heads _iT_k(z_i): UVΦscale/shift=Wbase(h)+∑k=1K[Wbasis(h)]k⋅T¯kUV_ _scale/shift=W_base(h)+ _k=1^K[W_basis(h)]_k· T_k (17) Discrete Cosine Transform (DCT) Kernel: The DCT kernel captures spectral dynamics across diffusion steps, it first forms a length-L signal u by combining frequency-damped spectral weights from h (with learnable decay exponent α via 1/(f+1)softplus(α)+ϵ1/(f+1)^softplus(α)+ε where ϵ∈(0,1)ε∈(0,1) is a constant added for numerical stability) with a fixed orthonormal Type-I DCT cosine basis ϕf\ _f\, where f indexes cosine frequency bins [15], then applies the same bounded structure as (3)–(4) with readouts from u, and for the B,CB,C base path, from h as well (≈ 7.1M parameters): u=∑f[Wspech]f(f+1)softplus(α)+ϵϕf.where ϵ=1×10−4u= _f [W_spech]_f(f+1)^softplus(α)+ε\, _f. $ε=1× 10^-4$ (18) Diagonal shifts and low-rank U,VU,V factors are Θshift _shift =WΘu, =W_ u, (19) UVΦscale/shift UV_ _scale/shift =WΦbaseh+WΦbasisu, =W_ _baseh+W_ _basisu, (20) 4 Experiments and Analysis 4.1 AQS Certificate and Diagnostics The Affine Quadratic Stability (AQS) certificate provides an analytic stability check for the bounded diagonal recurrence induced by the MaRK kernels [16]. Rather than solving a generic Semi-Definite Program (SDP) to search for a common Lyapunov function, we exploit the constructive parameterization chain from §2.1 to obtain an analytical certificate with P=IP=I. We load the 23-layer pretrained Hydra base model and extract each layer’s diagonal recurrence parameters AlogA_ . Since the maximum modulation bound βA _A is shared across the certified bounded variants, the bounded hyper-rectangle [Alog−βA,Alog+βA][A_ - _A,\;A_ + _A] is kernel-independent within this setting. For each layer ℓ and head i, we compute the worst-case discrete eigenvalue at the vertex δ=−βAδ=- _A, where A¯i A_i is closest to 11 (monotonicity of A¯i(δ) A_i(δ) in δ, see Proposition 2.1), using the minimum discretization step Δmin=0.001 _ =0.001. The energy function V(x)=‖x‖2V(x)=\|x\|^2 (i.e., P=IP=I) satisfies: A¯⊤IA¯−I=diag(A¯i2−1)≺0 A I\, A-I=diag( A_i^2-1) 0 (21) with per-layer stability margin εℓ=1−maxiA¯i,ℓ2>0 _ =1- _i A_i, ^2>0. All 23 layers certify with positive εℓ _ under the stated modulation and discretization assumptions, supporting stability of the modulated Hydra recurrence across admissible conditioning contexts. The per-layer spectral radii and stability margins are visualized in Figure 1. Figure 1: Stability margins (εℓ _ ) across all the layers of the pre-trained Hydra SSM at the extreme bounds of modulation by MaRK. 4.2 Synthetic LPV Identification This experiment evaluates the learnability and theoretical identifiability of the MaRK architecture. While the internal state matrices inhabit an unidentifiable manifold due to non-unique factorizations, the underlying input-output dynamics are estimable under the bounded LPV assumptions below. We evaluate recovery using the coordinate-invariant Markov Operator Error (MOE): ‖ℋtrue−ℋ^est‖F\|H_true- H_est\|_F (where ℋtrueH_true is the true Markov Parameter Sequence and ℋ^est H_est is the estimated Markov Parameter Sequence by MaRK). We train synthetic Hypernet, Chebyshev, and DCT adapter variants on N∈[50,1600]N∈[50,1600] trajectories with Gaussian output noise, keeping the backbone frozen to mirror our training regime; full synthetic protocol details are given in Appendix C. Proposition 4.1 (Operator Identifiability and Convergence). Let ℳM be an Affine Quadratically Stable LPV-SSM with strictly bounded representations, observed under standard excitation and mixing conditions. The precise statistical model, including the frozen-context identification protocol, bounded-density innovation assumptions, and operator-level persistent excitation, is specified in Appendix B.2. Then (i) the input-output operator defined by the frozen-time Markov sequence ℋ(t)=C′¯(t)A′k¯(t)B′¯(t)k=0∞H(c_t)=\ C (c_t) A \,k(c_t) B (c_t)\_k=0^∞ is identifiable up to a structural equivalence class. (i) Under squared loss, the empirical estimation error converges at the optimal parametric rate O(1/N)O(1/ N), where N is the total number of samples. Proof Sketch. We define model equivalence by equality of the input-output Markov sequence, which removes the non-identifiability of internal state coordinates while preserving the observable operator. The convergence argument then follows three standard steps. (i) Mixing: AQS gives geometric decay of the Markov operators; with persistently exciting conditioning and sub-Gaussian inputs/noise, as used in the synthetic recovery protocol, this places the induced temporal process in the exponentially β-mixing regime [21]. (i) Complexity: The tanh -bounded parameterization makes the operator class uniformly bounded and Lipschitz, so Talagrand’s contraction lemma yields Rademacher complexity scaling as O(1/N)O(1/ N) under the mixing relaxation of i.i.d. sampling [22]. (i) Generalization: Standard squared-loss bounds for dependent sequences transfer this capacity bound to the empirical operator-estimation error [23]. Appendix B.2 gives the full proof, including the structural equivalence relation and the cited statistical premises. ∎ (a) (b) (c) Figure 2: Markov Operator Error vs Dataset Size for synthetic Hypernet, Chebyshev, and DCT LPV recovery. Solid lines denote 15-seed means with normal-approximation 95% confidence intervals (shaded). Red dashed lines show an O(N−1/2)O(N^-1/2) reference normalized at the smallest dataset size. Analysis of Convergence. Figure 2 illustrates the log-log scaling of the MOE. All three synthetic kernels exhibit systematic operator-error decay as the number of trajectories increases, supporting the predicted O(1/N)O(1/ N) estimation trend in Proposition 4.1 under matched LPV recovery conditions. Fitting the endpoints over N=50N=50 to 16001600 gives approximate empirical rates of N−0.54N^-0.54 for Hypernet, N−0.66N^-0.66 for Chebyshev, and N−0.64N^-0.64 for DCT, consistent with parametric convergence up to kernel-dependent constants. Chebyshev achieves the lowest structured-basis error after the smallest sample regime, matching its smooth polynomial bias. DCT has a larger absolute constant, plausibly reflecting the additional bandlimited spectral restriction, but its monotone decay indicates that this restriction does not destroy operator-level recoverability. The experiment therefore supports the ability of MaRK kernel geometries to recover the coordinate-invariant Markov operator in a well-specified synthetic setting, while the theorem supplies a conditional statistical guarantee. 4.3 Dynamics Analysis To test whether MaRK conditions the temporal operator, we measure the Markov parameter sequence induced by each frozen diffusion timestep. For each adapter family, we evaluate the head-averaged Markov norm ‖hk(t)‖2\|h_k(t)\|_2 over lags k≤4096k≤ 4096 at t∈0,0.25,0.5,0.75,1.0t∈\0,0.25,0.5,0.75,1.0\, and compare it with the unmodulated Hydra backbone. Since ℋ(t)=C(t)Ak(t)B(t)k≥0H(t)=\C(t)A^k(t)B(t)\_k≥ 0 is the input-output memory profile of the SSM, separation between these curves provides evidence of timestep-dependent temporal kernels rather than merely timestep-dependent activations. The experiment therefore asks two questions: whether MaRK can reshape the Markov sequence as the diffusion process changes, and whether the resulting kernels remain dissipative at long lags. Figure 3: Head-resolved mean Markov parameter norm versus lag for Hypernet, Chebyshev, and DCT MaRK adapters at different diffusion timesteps t. The dashed curve denotes the frozen Hydra backbone without timestep-conditioned operator modulation. Analysis of Markov Dynamics. Figure 3 shows that all three adapter geometries preserve a stable decay envelope while producing distinct timestep-conditioned Markov sequences. The curves remain approximately geometric over thousands of lags, consistent with the AQS construction, but they do not collapse onto the frozen backbone: MaRK changes how much input-output energy is retained at each temporal scale. Hypernet produces the broadest upward shift relative to the backbone, especially at medium and long lags, suggesting an expressive increase in memory persistence. Chebyshev yields a tightly ordered family of smooth decay curves, matching the low-oscillation bias of polynomial modulation and indicating controlled, stable variation across diffusion time. DCT exhibits the most selective redistribution of energy: high-noise timesteps amplify early-lag response, while lower timesteps attenuate more rapidly and cross the backbone at intermediate lags, consistent with spectral shaping rather than uniform gain scaling. Together, these signatures support the claim that MaRK conditions the Markov parameter sequence, giving the diffusion timestep explicit influence over the SSM memory kernel while preserving long-horizon dissipation. Remark on Operator Adaptation. Figure 3 also illustrates all variants of MaRK possess a distinct signature compared to the frozen Hydra backbone. Since the Markov parameter sequence ℋkk=0∞\H_k\_k=0^∞ is the core characteristic of the SSM, it serves as a canonical operator for the system’s dynamics. This is because the sequence uniquely determines the transfer function and the impulse response of the model [18], remaining invariant under any similarity transformation of the state-space coordinates [24]. Consequently, the observed divergence in the norm signature across k-lags is consistent with a shift in the model’s impulse response. This shift suggests that MaRK does not merely perform a linear re-projection of features, but also adapts the underlying dynamical operators (specifically the spectral properties and memory decay) to align with the temporal dependencies of the target task. This supports viewing our approach as a principled mechanism for operator-level fine-tuning. 4.4 Ablations We ablate the three MaRK adapter geometries under the same frozen Hydra backbone, diffusion training objective, and CART-weighted validation protocol. The goal is to isolate the effect of the operator-modulation map ψ(⋅)ψ(·): direct Hypernet projections, smooth Chebyshev polynomial bases, or bandlimited DCT bases. Table 1 reports ℒCARTL_CART, where lower values indicate better reconstruction of masked tokens under the same context-adaptive weighting used for training. We include trainable adapter size to separate geometry effects from parameter count, and interpret the results as an internal MaRK design ablation rather than as a full scaling comparison against separately trained diffusion language models. Table 1: CART losses (ℒCART↓L_CART ) across 10 seeds (mean ± 95% CI, 1.96×1.96×SEM) and trainable adapter sizes across completed validation datasets. Bold = Best. Underline = Second Best. Kernel Params PTB WikiText Lambada AG News Arxiv Avg. Hypernet 6.3M 3.39 ± 0.31 3.73 ± 0.23 3.70 ± 0.17 3.91 ± 0.20 4.12 ± 0.10 3.77 Chebyshev 11.0M 2.42 ± 0.26 2.55 ± 0.16 2.40 ± 0.10 2.85 ± 0.14 2.55 ± 0.05 2.55 DCT 7.1M 2.48 ± 0.27 2.55 ± 0.14 2.47 ± 0.11 2.86 ± 0.14 2.61 ± 0.07 2.59 Both structured bases improve on the direct Hypernet projection on every completed validation dataset. Chebyshev gives the lowest average loss (2.55) and leads DCT on four of the five datasets, tying it on WikiText; DCT is 0.04 nats behind on the average (2.59). The cause is structural rather than statistical. The timestep is a continuous scalar and adjacent timesteps define nearly identical denoising tasks, so the modulation map should vary smoothly, but the Hypernet MLP consumes each timestep as an independent input and must learn that smoothness from data. It can therefore fit the training timesteps accurately while varying freely between them, spending capacity on high-frequency structure the diffusion process does not need. Chebyshev and DCT encode the smoothness in the basis itself: five global polynomials, or eight frequencies with a learned spectral decay that attenuates the high end. Projecting the modulation map onto a compact basis functions set acts as a strong structural regularizer, constraining the functional space without compromising expressive capacity. Table 1 shows the effect at equal budget: Hypernet and DCT carry 6.3M and 7.1M adapter parameters and still differ by 1.18 nats, while the further 3.9M that separate DCT from Chebyshev (11.0M) buys an additional 0.04 nats. Hypernet does not overfit; it lacks the prior. Chebyshev then separates from DCT only on Arxiv (2.55 vs. 2.61), Lambada (2.40 vs. 2.47) and PTB (2.42 vs. 2.48), with WikiText and AG News inside seed noise, and DCT is the cheaper of the two at inference, so Chebyshev is the accuracy ceiling and DCT the better parameter and latency trade. Table 2: Parameter ablation: CART validation loss (ℒCART↓L_CART ) on the WikiText dataset across 10 seeds (mean ± 95% CI, 1.96×1.96×SEM). "Full" modulates all five parameters; "None" modulates none (frozen Hydra). Bold = Best per kernel. Kernel Full All Except A A Only Δ Only B,C Only D Only All Except Δ None Hypernet 3.73 ± 0.26 3.91 ± 0.27 4.07 ± 0.29 4.04 ± 0.29 4.08 ± 0.29 3.78 ± 0.26 3.74 ± 0.26 3.99 ± 0.29 Chebyshev 2.55 ± 0.18 2.78 ± 0.19 3.44 ± 0.22 3.92 ± 0.26 3.06 ± 0.21 3.90 ± 0.26 2.56 ± 0.18 4.00 ± 0.27 DCT 2.55 ± 0.16 2.80 ± 0.16 3.24 ± 0.17 3.99 ± 0.24 3.05 ± 0.17 3.90 ± 0.24 2.59 ± 0.16 4.00 ± 0.24 A full leave-one-out parameter ablation across all five state-space model (SSM) parameters (A,B,C,D,ΔA,B,C,D, ) isolated the performance contributions of each operator-modulation component. Modulating matrix A proves necessary across all settings; freezing A (All Except A) increases validation loss by +0.18 to +0.25 points. This performance gap confirms that MaRK cannot be reduced to standard Mamba-style selection (B,C,ΔB,C, ), as the B,CB,C-Only variant leaves a substantial loss degradation of +0.34 to +0.51 points. This result aligns with the functional distinction where input gating controls instantaneous throughput while A governs persistent memory horizon [12, 1]. Furthermore, modulating Δ yields marginal independent gains due to the redundancy in spectral transformation paths: the operator dynamics can be shifted directly via A or indirectly via Δ through zero-order hold discretization A¯=exp(Acont⋅Δ) A= (A_cont· ). The direct spectral path dominates empirically, as A-Only outperforms Δ -Only by up to 0.75 points, and All Except Δ achieves the second-best score across every kernel. Nevertheless, single-parameter configurations trail the Full baseline by up to +0.50 points, demonstrating that optimal performance requires joint parameter modulation. Table 3: Inference performance and relative latency overhead over Frozen Hydra at batch size 8 on an NVIDIA A100-SXM4-80GB with the same frozen 23-layer backbone. Each variant: 15 seeds, 100 forward passes per seed, means ± 95% CI (t-distribution, 14df). Bold = MaRK (Ours) Variant Latency (ms) seq=4096 Memory (MiB) seq=4096 Total Params (M) Additional Params (M) Overhead seq=512 Overhead seq=4096 Frozen Hydra 267.43 ± 0.06 4446.58 ± 0.06 111.68 0.00 - - Input Injection 270.86 ± 0.03 4455.25 ± 0.06 113.96 2.28 +3.6% +1.3% AdaLN Zero 276.84 ± 0.02 4474.52 ± 0.60 118.51 6.84 +7.6% +3.5% Hypernet 304.75 ± 0.03 4470.62 ± 0.06 117.95 6.27 +65.6% +14.0% Chebyshev 309.74 ± 0.10 4488.61 ± 0.06 122.65 10.97 +104.0% +15.8% DCT 307.28 ± 0.05 4473.93 ± 0.06 118.75 7.07 +88.7% +14.9% At sequence length 4096, MaRK incurs a modest 14% to 16% latency increase while maintaining memory overhead below 1% relative to the frozen Hydra backbone. The per-layer conditioning computation introduces a constant temporal additive cost (∼ 37 ms to 42 ms) that remains independent of sequence length. Consequently, while conditioning overhead dominates leading to 65% to 104% additional latency at length 512, this cost amortizes rapidly as sequence length scales to 4096 and beyond. Because MaRK scales O(1)O(1) with sequence length and O(N)O(N) with diffusion sampling steps, its asymptotic runtime trajectory aligns with standard AdaLN Zero and Input-Injection baselines. 5 Conclusion We introduced Dynamic Operator Modulation, an operator-level conditioning perspective that shifts SSM adaptation from external signal injection toward direct reconfiguration of the recurrent operator. By mapping context to bounded, low-rank modulations of the frozen backbone, MaRK induces a context-indexed family of stable input-output dynamics. This operator-level perspective is our central contribution: adapting an SSM for iterative tasks like diffusion can be understood as a controlled deformation of the model’s underlying Markov parameters. Our framework provides theoretical support (through Affine Quadratic Stability certificates and Rademacher complexity bounds) and empirical evidence across multiple adapter geometries. The observed shift in Markov signatures suggests that parameter-efficient fine-tuning can arise from the structure of this modulation, rather than only from an additive patch. Overall, these results suggest that the temporal logic of an SSM can be made context-aware in the studied setting while preserving stability and recoverability under the stated assumptions, opening new avenues for principled adaptation in non-stationary sequence modeling. Systematic empirical comparison against established conditioning frameworks on downstream benchmarks, and how MaRK holds under scaling laws remain important future directions; we discuss these limitations further in Appendix A.5. 6 Acknowledgements This project was made possible by the generous support provided by the NVIDIA AI Technology Center (NVAITC). Additionally, we gratefully acknowledge Dr. Yuhao Wang from the Department of Data Science at City University of Hong Kong for his assistance during the final stages of this project. References [1] A. Gu and T. Dao (2024) Mamba: linear-time sequence modeling with selective state spaces. In First Conference on Language Modeling, External Links: Link Cited by: §1, §4.4. [2] T. Dao and A. Gu (2024) Transformers are ssms: generalized models and efficient algorithms through structured state space duality. In Forty-first International Conference on Machine Learning, ICML 2024, Vienna, Austria, July 21-27, 2024, R. Salakhutdinov, Z. Kolter, K. A. Heller, A. Weller, N. Oliver, J. Scarlett, and F. Berkenkamp (Eds.), Proceedings of Machine Learning Research, p. 10041–10071. External Links: Link Cited by: §1. [3] A. Lahoti, K. Y. Li, B. Chen, C. Wang, A. Bick, J. Z. Kolter, T. Dao, and A. Gu (2026) Mamba-3: improved sequence modeling using state space principles. CoRR abs/2603.15569. External Links: Link, Document, 2603.15569 Cited by: §1. [4] J. Ho, A. Jain, and P. Abbeel (2020) Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual, H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin (Eds.), External Links: Link Cited by: §1, §2. [5] J. Austin, D. D. Johnson, J. Ho, D. Tarlow, and R. van den Berg (2021) Structured denoising diffusion models in discrete state-spaces. In Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual, M. Ranzato, A. Beygelzimer, Y. N. Dauphin, P. Liang, and J. W. Vaughan (Eds.), p. 17981–17993. External Links: Link Cited by: §1, §3.2. [6] X. L. Li, J. Thickstun, I. Gulrajani, P. Liang, and T. B. Hashimoto (2022) Diffusion-lm improves controllable text generation. In Advances in Neural Information Processing Systems 35: Annual Conference on Neural Information Processing Systems 2022, NeurIPS 2022, New Orleans, LA, USA, November 28 - December 9, 2022, S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (Eds.), External Links: Link Cited by: §1. [7] S. Nie, F. Zhu, Z. You, X. Zhang, J. Ou, J. Hu, J. ZHOU, Y. Lin, J. Wen, and C. Li (2025) Large language diffusion models. In ICLR 2025 Workshop on Deep Generative Model in Machine Learning: Theory, Principle and Efficacy, External Links: Link Cited by: §1, §3.2, §3.2, §3.2. [8] J. Ye, Z. Xie, L. Zheng, J. Gao, Z. Wu, X. Jiang, Z. Li, and L. Kong (2025) Dream 7b: diffusion large language models. CoRR abs/2508.15487. External Links: Link, Document, 2508.15487 Cited by: §A.3, §A.3, §1, §1, §3.2, §3.2, §3.2, §3.3. [9] W. Peebles and S. Xie (2023) Scalable diffusion models with transformers. In IEEE/CVF International Conference on Computer Vision, ICCV 2023, Paris, France, October 1-6, 2023, p. 4172–4182. External Links: Link, Document Cited by: §1, §2. [10] J. N. Yan, J. Gu, and A. M. Rush (2024) Diffusion models without attention. In IEEE/CVF Conference on Computer Vision and Pattern Recognition, CVPR 2024, Seattle, WA, USA, June 16-22, 2024, p. 8239–8249. External Links: Link, Document Cited by: §1. [11] H. Phung, Q. Dao, T. T. Dao, H. Phan, D. N. Metaxas, and A. T. Tran (2024) DiMSUM: diffusion mamba - a scalable and unified spatial-frequency method for image generation. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, External Links: Link Cited by: §1. [12] A. Gu, K. Goel, and C. Ré (2022) Efficiently modeling long sequences with structured state spaces. In The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022, External Links: Link Cited by: §1, §2, §4.4. [13] S. Hwang, A. S. Lahoti, R. Puduppully, T. Dao, and A. Gu (2024) Hydra: bidirectional state space models through generalized matrix mixers. In Advances in Neural Information Processing Systems 38: Annual Conference on Neural Information Processing Systems 2024, NeurIPS 2024, Vancouver, BC, Canada, December 10 - 15, 2024, A. Globersons, L. Mackey, D. Belgrave, A. Fan, U. Paquet, J. M. Tomczak, and C. Zhang (Eds.), External Links: Link Cited by: §1, §1, §3.1, §3.1, §3.1. [14] E. J. Hu, Y. Shen, P. Wallis, Z. Allen-Zhu, Y. Li, S. Wang, L. Wang, and W. Chen (2022) LoRA: low-rank adaptation of large language models. In The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022, External Links: Link Cited by: §1. [15] N. Ahmed, T. Natarajan, and K. R. Rao (1974) Discrete cosine transform. IEEE Transactions on Computers 100 (1), p. 90–93. Cited by: §1, §3.4. [16] P. B. Cox (2018) Towards efficient identification of linear parameter-varying state-space models. Ph.D. Thesis, Eindhoven University of Technology, Eindhoven, The Netherlands. Note: Proefschrift External Links: ISBN 978-90-386-4456-1, Link Cited by: 1st item, §1, §2, §4.1. [17] J. S. Shamma (1988) Analysis and design of gain scheduled control systems. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA. Note: The LPV-SSM equation is introduced on page 42. External Links: Link Cited by: §2. [18] T. Kailath (1980) Linear systems. Information and System Sciences Series, Prentice-Hall. Note: Markov parameter sequence Hk=CAk−1BH_k=CA^k-1B (for k≥1k≥ 1) defined on p. 70, 113 and 354. External Links: ISBN 9780135369616, LCCN 79014928, Link Cited by: §2, §4.3. [19] L. Zhu, B. Liao, Q. Zhang, X. Wang, W. Liu, and X. Wang (2024) Vision mamba: efficient visual representation learning with bidirectional state space model. In Forty-first International Conference on Machine Learning, ICML 2024, Vienna, Austria, July 21-27, 2024, R. Salakhutdinov, Z. Kolter, K. A. Heller, A. Weller, N. Oliver, J. Scarlett, and F. Berkenkamp (Eds.), Proceedings of Machine Learning Research, p. 62429–62442. External Links: Link Cited by: §3.1. [20] W. Chen, L. Niu, Z. Lu, F. Meng, and J. Zhou (2024) MaskMamba: A hybrid mamba-transformer model for masked image generation. CoRR abs/2409.19937. External Links: Link, Document, 2409.19937 Cited by: §3.1. [21] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan (1994) Linear matrix inequalities in system and control theory. Studies in Applied Mathematics, Vol. 15, Society for Industrial and Applied Mathematics (SIAM). External Links: ISBN 0-89871-334-X Cited by: §4.2. [22] M. Ledoux and M. Talagrand (1991) Probability in banach spaces: isoperimetry and processes. Springer-Verlag, Berlin, New York. External Links: ISBN 978-3-642-20212-4, Document Cited by: §B.2, §4.2. [23] P. L. Bartlett and S. Mendelson (2002) Rademacher and Gaussian complexities: risk bounds and structural results. Journal of Machine Learning Research 3 (Nov), p. 463–482. Cited by: §B.2, §4.2. [24] C. Chen (1998) Linear system theory and design. 3rd edition, The Oxford Series in Electrical and Computer Engineering, Oxford University Press, New York, NY. External Links: ISBN 9780195117776 Cited by: §4.3. [25] Y. Schiff, S. S. Sahoo, H. Phung, G. Wang, S. Boshar, H. Dalla-torre, B. P. de Almeida, A. M. Rush, T. PIERROT, and V. Kuleshov (2025) Simple guidance mechanisms for discrete diffusion models. In The Thirteenth International Conference on Learning Representations, External Links: Link Cited by: §A.3. [26] D. Ha, A. M. Dai, and Q. V. Le (2017) HyperNetworks. In 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings, External Links: Link Cited by: item (i), item (i), §A.4. [27] E. Haber and L. Ruthotto (2017) Stable architectures for deep neural networks. Inverse Problems 34 (1), p. 014004. External Links: ISSN 1361-6420, Link, Document Cited by: item (i). [28] T. Galanti and L. Wolf (2020) On the modularity of hypernetworks. In Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual, H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin (Eds.), External Links: Link Cited by: item (i). [29] T. D. Pham and L. T. Tran (1985) Some mixing properties of time series models. Stochastic Processes and their Applications 19 (2), p. 297–303. Cited by: §B.2, Lemma B.3. [30] B. Yu (1994) Rates of convergence for empirical processes of stationary mixing sequences. The Annals of Probability 22 (1), p. 94–116. External Links: Document, Link Cited by: §B.2. [31] M. Mohri and A. Rostamizadeh (2008) Rademacher complexity bounds for non-i.i.d. processes. In Advances in Neural Information Processing Systems 21, Proceedings of the Twenty-Second Annual Conference on Neural Information Processing Systems, Vancouver, British Columbia, Canada, December 8-11, 2008, D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou (Eds.), p. 1097–1104. External Links: Link Cited by: §B.2, §B.2. [32] V. Kuznetsov and M. Mohri (2017) Generalization bounds for non-stationary mixing processes. Machine Learning 106 (1), p. 93–117. Cited by: §B.2. Appendix A Discussions A.1 Notation Summary Table 4: Compact symbol guide for the main text. Symbol Meaning t,τt,τ diffusion timestep and internal continuous time tc_t context/timestep embedding used as the LPV scheduling variable Q,Q′Q,Q base and modulated SSM parameters, Q∈A,B,C,D,ΔQ∈\A,B,C,D, \ A(t),B(t)A(c_t),B(c_t) LPV-SSM matrix-valued functions for recurrence and read-in parameters C(t),D(t)C(c_t),D(c_t) LPV-SSM matrix-valued functions for read-out and skip parameters x(τ),u(τ),y(τ)x(τ),u(τ),y(τ) LPV-SSM vector-valued functions for state, input and output A,B,C,D,ΔA,B,C,D, SSM recurrence, read-in, read-out, skip, and discretization parameters Alog′,Δinput,Δbias′A_ , _input, _bias modulated log-recurrence, discretization input, and modulated discretization bias Θ,Φ , diagonal parameters A,Δ,D\A, ,D\ and dense projections B,C\B,C\ ψscale,ψshift _scale, _shift adapter maps that emit bounded scales and shifts αΦ,βΦ,βΘ _ , _ , _ modulation magnitudes for scales and shifts P,I,εP,I, Lyapunov matrix, identity matrix, and AQS stability margin A¯,B¯ A, B ZOH-discretized recurrence and read-in parameters ℋ(t),hkH(c_t),h_k context-indexed Markov parameter sequence and its kkth lag term ,ℳ,θV,M,θ vocabulary, absorbing mask token, and trainable model parameters h,W⋅,z,Tk,ϕfh,W_·,z,T_k, _f adapter hidden state, learned readouts, Chebyshev coordinate/polynomial, and DCT basis w(t),wi(t,t)w(t),w_i(t,x_t) global diffusion weight and CART token-level weight U,V,rU,V,r low-rank factors and rank used for dense-parameter modulation A.2 The MaRK Kernel Classes We propose three principal variants of the MaRK kernels: Hypernet, Chebyshev, and DCT. Each kernel structurally modifies the recurrent dynamics (A,ΔA, ) and input/output projections (B,CB,C) of the base State Space Model. In the certified bounded setting, shifts are modulated by bounded activation functions (hyperbolic tangent gating) to support AQS and prevent unbounded activation growth. The three variants should be read as different parameterizations of the same bounded operator-modulation map rather than as unrelated architectures. Hypernet places the weakest structural prior on the context-to-operator map and is therefore useful as a flexible baseline. Chebyshev and DCT instead restrict the modulation to structured bases: Chebyshev favors globally smooth variation over the conditioning coordinate, while DCT favors bandlimited spectral variation with explicit attenuation of higher-frequency components. This distinction is important because the empirical comparison in Table 1 tests not only whether MaRK can modulate the SSM operator, but also which inductive bias makes that modulation reliable under a frozen backbone and finite adaptation budget. • Expressivity. Hypernet can represent the least constrained modulation among the three variants, but this flexibility may allocate capacity to changes that are not necessary for the task. Chebyshev and DCT reduce this search space by forcing modulations to pass through low-dimensional basis expansions. • Optimization and parameter efficiency. The basis kernels trade some pointwise flexibility for smoother parameter trajectories. This can make the adapter easier to regularize and interpret, although the actual compute and parameter cost still depends on basis degree, frequency budget, and low-rank factor sizes. • Interpretability. Chebyshev coefficients can be interpreted as polynomial components of the operator path, while DCT coefficients separate low- and high-frequency modulation content. Hypernet is less directly interpretable, since the MLP output mixes these effects before the bounded shifts are applied. A.2.1 Hypernet Adapter The Hypernet acts as the fundamental MaRK baseline, employing standard linear projections directly from the hidden conditioning variable to the factorized shift components. This variant serves as the most expressive baseline because the MLP can learn a comparatively unconstrained mapping from the conditioning vector to the modulations. Its main role in the study is therefore diagnostic: it asks whether direct learned projections are sufficient, before imposing the smoother basis restrictions used by Chebyshev or DCT. Algorithm 1 Hypernet MaRK Kernel 1: Input: Context vector tc_t, Base parameters Alog,B,C,Δ,DA_ ,B,C, ,D 2: Hyperparameters: Factorization constraints α,βα,β 3: h←MLP(t)h (c_t) ⊳ Extract hidden representation 4: Ashift,Δshift,Dshift←LinearA,Δ,D(h)A_shift, _shift,D_shift _A, ,D(h) 5: Bscale,Bshift←LowRankProject(h)B_scale,B_shift (h) 6: Cscale,Cshift←LowRankProject(h)C_scale,C_shift (h) 7: Modulation: 8: Alog′←Alog+βAtanh(Ashift)A_ ← A_ + _A (A_shift) 9: Δ′←Δ+βΔtanh(Δshift) ← + _ ( _shift) 10: B′←B⊙(1+αBtanh(Bscale))+βBtanh(Bshift)B ← B (1+ _B (B_scale))+ _B (B_shift) 11: C′←C⊙(1+αCtanh(Cscale))+βCtanh(Cshift)C ← C (1+ _C (C_scale))+ _C (C_shift) 12: Return modulated parameters A.2.2 Chebyshev Polynomial Adapter The Chebyshev polynomial adapter imposes structural smoothness by mapping conditional modulations onto orthogonal polynomial bases. Rather than predicting unbounded shifts directly, the hidden dimension governs a latent evaluation coordinate z∈[−1,1]z∈[-1,1] and generates series coefficients. This creates a low-degree global approximation to the operator path, which is appropriate when the desired recurrent dynamics vary smoothly with context. The same global support also means that changing one coefficient can affect the modulation over the full coordinate range, so the degree K controls a bias-variance tradeoff rather than merely increasing capacity. Algorithm 2 Chebyshev MaRK Kernel 1: Input: Context vector tc_t, Base parameters Alog,B,C,Δ,DA_ ,B,C, ,D 2: Hyperparameters: Degree K, Scaling constraints α,βα,β 3: h←MLP(t)h (c_t) 4: z←Linearz(h)z _z(h) ⊳ Project hidden state to latent basis coordinate 5: T(z)←[T0(z),T1(z),…,TK−1(z)]T(z)←[T_0(z),T_1(z),…,T_K-1(z)] ⊳ Evaluate Chebyshev polynomials up to degree K 6: θcoeffs←Linearθ(h)for θ∈A,Δ,D _coeffs _θ(h) θ∈\A, ,D\ 7: θshift←∑k=0K−1θcoeffs(k)⋅Tk(z)for θ∈A,Δ,D _shift← _k=0^K-1 _coeffs^(k)· T_k(z) θ∈\A, ,D\ 8: T¯←1nheads∑iT(zi) T← 1n_heads _iT(z_i) ⊳ Calculate shared basis mean for broad matrix shifts 9: Xscale←UVXscale_base(h)+∑k=0K−1UVXscale_basis(h)(k)⋅T¯kfor X∈B,CX_scale _X_scale\_base(h)+ _k=0^K-1UV_X_scale\_basis(h)^(k)· T_k X∈\B,C\ 10: Xshift←UVXshift_base(h)+∑k=0K−1UVXshift_basis(h)(k)⋅T¯kfor X∈B,CX_shift _X_shift\_base(h)+ _k=0^K-1UV_X_shift\_basis(h)^(k)· T_k X∈\B,C\ 11: Modulation: Evaluate modulated representations using base parameters and Xscale/shiftX_scale/shift via bounded tanh saturation (Step 8-11 in Algorithm 1) 12: Return modulated parameters A.2.3 Discrete Cosine Transform (DCT) Kernel The DCT kernel imposes a spectral inductive bias by casting modulations as bandlimited real-valued basis transformations. Rather than assuming that high-frequency modulation is always harmful, the parameterization makes high-frequency content explicit and attenuates it through a learned monotone decay. Drawing inspiration from the Inverse Real Fast Fourier Transform (IRFFT) without requiring unconstrained complex arithmetic, we formulate a deterministic real-valued basis grid. Contextual representations evaluate spectral components that decay across higher frequencies before being mapped back into parameter shifts. This makes DCT most natural when the useful operator variation is expected to be coarse-to-fine or spectrally organized; if the true modulation requires sharp localized changes, the frequency budget and decay prior can become limiting factors. Algorithm 3 DCT MaRK Kernel 1: Input: Context vector tc_t, Base parameters Alog,B,C,Δ,DA_ ,B,C, ,D 2: Hyperparameters: NumFreqs nfn_f, Basis grid timepoints L, Decay bound αdecay _decay 3: h←MLP(t)h (c_t) 4: wdecay(f)←1(f+1)softplus(α)+ϵw_decay(f)← 1(f+1)^softplus(α)+ε ⊳ Evaluate learned frequency attenuation 5: Specweights←Linearspec(h)⊙wdecaySpec_weights _spec(h) w_decay 6: Timevec←Specweights×ℬDCT−IITime_vec _weights×B_DCT-I ⊳ Project bandlimited spectral amplitudes into temporal structure 7: θshift←Linearθ(Timevec)for θ∈A,Δ,D _shift _θ(Time_vec) θ∈\A, ,D\ 8: Xscale←UVXscale_base(h)+UVXscale_basis(Timevec)for X∈B,CX_scale _X_scale\_base(h)+UV_X_scale\_basis(Time_vec) X∈\B,C\ 9: Xshift←UVXshift_base(h)+UVXshift_basis(Timevec)for X∈B,CX_shift _X_shift\_base(h)+UV_X_shift\_basis(Time_vec) X∈\B,C\ 10: Modulation: Evaluate modulated representations using base parameters and Xscale/shiftX_scale/shift via bounded tanh saturation (Step 8-11 in Algorithm 1) 11: Return modulated parameters Overall, the kernel choice determines how much structure is imposed before the bounded AQS-compatible shifts are applied. Hypernet may overfit or produce rapidly varying modulations when the available data do not justify its flexibility. Chebyshev may oversmooth localized changes because its polynomial basis has global support. DCT may underfit non-bandlimited operator paths if too few frequencies are retained or if the learned decay is too strong. These are not separate stability assumptions, since all variants still use bounded shifts, but they are practical modeling assumptions that affect how efficiently MaRK uses its trainable degrees of freedom. A.3 CART Objective To rigorously adapt the training criteria against distributional mode collapse in heavily masked regimes, we replace standard Categorical Cross-Entropy (CCE) masked language modeling loss with the Context-Adaptive Token-Level Noise Rescheduling (CART) objective. In conventional sequence modeling, heavy masking at early diffusion steps often degenerates the output tensor into identical redundant tokens predicting the global marginal distribution rather than contextualized completions. To prevent mode collapse, the objective introduces differentiated spatial penalty terms, adapted from Dream [8]. As a consequence, when evaluating with a global context such as MDLM ELBO [25], the perplexities may seem inflated. We provide the full mathematical specification of the Context-Adaptive Token-level Noise Rescheduling (CART) objective described in §3.3, including the token-wise weighting function and its substitution into the masked diffusion loss. Let mi=[xt,i=ℳ]m_i=1[x_t,i=M] indicate which tokens are masked at timestep t. Following Dream [8], CART replaces the single sequence-level weight with a per-token weight built from a symmetric geometric kernel over token distance. For a masked token i, wi(t,t)=12∑j=1L(1−mj)Geo(p,|i−j|−1),Geo(p,d)=p(1−p)d,w_i(t,x_t)= 12 _j=1^L(1-m_j)\;Geo\! (p,\,|i-j|-1 ), (p,d)=p\,(1-p)^d, (22) where the sum runs over the visible (unmasked) tokens and p∈(0,1]p∈(0,1] controls the kernel sharpness. A small p spreads each clean token’s influence almost uniformly across the sequence, while a large p concentrates it on nearby masked tokens. Because the weight is the sum of one geometric kernel per visible token, a masked position surrounded by many nearby visible tokens receives a larger weight; this is the mixture of geometric kernels over clean anchors used in Dream. We use p=0.45p=0.45 (Appendix C.2). Substituting wi(t,t)w_i(t,x_t) into Eq. (13) yields the CART-weighted training objective: ℒCART(θ)=0,t,t[∑i=1Lmiwi(t,t)(−logpθ(x0,i∣t,t))].L_CART(θ)=E_x_0,\,t,\,x_t [ _i=1^Lm_i\,w_i(t,x_t)\; (- p_θ(x_0,i _t,t) ) ]. (23) A.4 Scope of Contribution Before discussing limitations, we clarify the scope of what this paper establishes. Our contribution is a theoretically grounded formulation of Dynamic Operator Modulation for SSMs: instead of injecting conditioning information only through inputs or activation scaling, MaRK maps a context vector into bounded changes of the recurrent operator itself. The comparison to Adaptive Layer Normalization (AdaLN) and input-stream injection in Section 1 is a structural motivation, not a claim that MaRK dominates those mechanisms in all empirical settings. Scope of the conditioning signal. The operator-conditioning claim is scoped to a continuous, scalar, ordered scheduling variable, which in this work is the diffusion timestep t∈[0,1]t∈[0,1]. This setting matters because adjacent timesteps index nearly identical denoising tasks, so the context-to-operator map should vary smoothly with the coordinate. Chebyshev and DCT encode that smoothness structurally, which is the mechanism behind their advantage over the Hypernet adapter (Section 4.4). The MaRK architecture accepts any embedding as tc_t, but the smooth-operator prior is specific to a continuous and ordered coordinate, so we do not claim transfer to autoregressive next-token modeling, which has no denoising coordinate to schedule against, or to high-dimensional prompt conditioning, where the signal is discrete and carries no proximity metric. Those regimes are left to future work. The central technical challenge is that Dynamic Operator Modulation is a meta-function learning problem: a conditioning vector tc_t must determine continuous-time system parameters (A,B,C,D,Δ)(A,B,C,D, ) that control long-range temporal behavior. This resembles the setting studied by hypernetworks [26], where a side network produces parameters for a primary model. In recurrent or state-space systems, such context-to-operator maps raise three practical risks: (i) Training instability and divergence. If generated recurrence parameters are unconstrained, small context-dependent errors can move eigenvalues toward unstable regions, producing unstable forward dynamics or poorly conditioned gradients [26, 27]. (i) Identifiability collapse. A low-dimensional conditioning signal can map to many internally different state-space parameterizations that represent the same input-output behavior, and poorly structured maps can become nearly insensitive to tc_t [28]. (i) Optimization burden. Emitting large dense parameter tensors from a small context vector can introduce high-variance gradients and make the context-to-operator map difficult to train reliably [26]. MaRK addresses these risks through a constrained architecture rather than by generating unrestricted recurrent weights. • Stability is controlled by bounded log-space modulations and positive discretization, yielding the constructive AQS certificate in Section 4.1 and Appendix B under the stated parameterization [16]. • Operator recovery is studied at the level of the coordinate-invariant Markov parameter sequence. The synthetic LPV experiments support recovery of this input-output operator under matched assumptions, while the theory states the corresponding assumptions explicitly. • Trainability is encouraged by low-dimensional outputs: the adapters predict scalar shifts and low-rank scales rather than full dense matrices, using 6.3–11M trainable parameters on top of a frozen 111M-parameter backbone. Thus, the contribution of this work is the formulation and initial validation of an operator-level conditioning paradigm for SSMs. The parameter ablation in Table 2, the inference-overhead comparison in Table 3, and the Markov-operator diagnostics in Section 4.3 provide mechanistic and systems-level evidence that the conditioning signal reaches the recurrent operator rather than only external activations. None of these is a matched-budget quality comparison against external conditioning methods, and Appendix A.5 states what remains open. A.5 Limitations and Future Directions This work establishes MaRK as a stable and identifiable operator-level conditioning framework for SSMs, but its empirical scope is intentionally narrower than the full space of conditioning methods. Our experiments validate internal MaRK design choices (kernel geometry in Table 1, SSM-parameter coverage in Table 2), compare systems cost against external conditioning methods (Table 3), and diagnose changes in the Markov parameter sequence (Figure 3). The comparison against AdaLN and input-stream injection is limited to inference latency, memory, and parameter count. We do not provide a controlled head-to-head baseline against those methods under matched training compute, objective, and frozen-backbone constraints, so we make no claim that operator modulation outperforms them on reconstruction quality. The present results also do not establish scaling laws for Dynamic Operator Modulation. All reported language-model experiments use a 111M-parameter Hydra-BERT backbone with 6.3–11M trainable adapter parameters. This setting is sufficient to test whether bounded operator modulation is feasible, stable, and trainable, but it does not determine how MaRK behaves as model depth, hidden width, adapter rank, sequence length, training tokens, or pretrained SSM family vary. In particular, the Markov-dynamics diagnostics in Figure 3 suggest that MaRK can reshape long-lag memory profiles, but they should be interpreted as mechanistic evidence rather than as a completed long-context scaling study. Our evaluation is likewise centered on the claims of this paper: AQS certificates, synthetic LPV recovery, Markov-operator diagnostics, and CART-weighted validation loss. These measurements support stability, recoverability under matched assumptions, and timestep-conditioned operator adaptation. We now report inference systems metrics (Table 3), but the evaluation does not cover generation quality, downstream task accuracy, or robustness across long-context language benchmarks. Future work should test MaRK on controlled long-range suites such as LRA and language long-context evaluations, and should report both quality and systems metrics there. The theoretical guarantees depend on assumptions that may be violated in less structured settings. The identifiability and generalization arguments assume a well-specified bounded LPV-SSM class, stable Markov operators, sufficiently exciting inputs, sub-Gaussian noise, and mixing conditions. The written stability proof relies on the bounded MaRK parameterization and constructive ZOH chain, while the statistical-learning proof composes standard concentration, mixing, and Rademacher-complexity results cited from the literature. These assumptions are appropriate for the synthetic recovery protocol and for the stability claims made here, but they should not be read as unconditional guarantees for arbitrary pretrained sequence models or arbitrary conditioning distributions. Finally, the bounded parameterization that makes MaRK stable may also constrain expressivity. Performance can plausibly depend on the modulation bounds, the low-rank dimension, the Chebyshev degree, the DCT frequency count, and the diffusion weighting. A useful next stage is a systematic robustness study over these hyperparameters. The conditioning regime itself is also open: within this paper the scheduling variable is the diffusion timestep, and whether operator-level conditioning helps when the conditioning signal is not a continuous and ordered coordinate remains future work (Appendix A.4). Appendix B Proofs This section gives the written theoretical arguments for the strictly bounded MaRK kernels. We first state the AQS stability argument for bounded affine parameter variations, then prove the coordinate-invariant operator identifiability and statistical convergence guarantee for Proposition 4.1. B.1 Proof of Affine Quadratic Stability (AQS) Let the base matrix A∈ℝn×nA ^n× n and perturbed system A(v)=A+ΔA(v)A(v)=A+ A(v), where parameter variations v lie in a hyper-rectangle Ψ . Theorem. The LPV system is Affine Quadratically Stable if there exists a matrix P≻0P 0 and scalar ε>0 >0 such that for all vertices vk∈Ψv_k∈ , A(vk)⊤PA(vk)−P⪯−εIA(v_k) PA(v_k)-P - I. Proof. Because the mapping t↦A(v)c_t A(v) is modulated via strictly bounded differentiable activations (tanh ), the range of possible evaluations falls strictly within the convex hull defined by vertices vkv_k. This vertex argument applies to affine dependence on the bounded modulation variables; in the diagonal Hydra ZOH instantiation, the certificate is obtained directly with P=IP=I because A¯i=exp(−exp(Alog,i′)Δi)∈(0,1) A_i= (- (A _ ,i) _i)∈(0,1). For any affine interpolation A(λ)=∑kλkAkA(λ)= _k _kA_k, matrix convexity of A↦A⊤PA A PA for P≻0P 0 gives A(λ)⊤PA(λ)⪯∑kλkAk⊤PAk.A(λ) PA(λ) _k _kA_k PA_k. Thus, if the strictly negative definite LMI Ak⊤PAk−P⪯−εIA_k PA_k-P - I is satisfied uniformly at all finite vertices, every interpolated system inherits the same Lyapunov decrease. Hence V(x)=x⊤PxV(x)=x Px decays exponentially. ∎ B.2 Proof of Operator Identifiability (Proposition 4.1) We prove the result for the coordinate-invariant input-output operator. The direct term Dθ(c)D_θ(c) is not part of the state-coordinate ambiguity, but it is part of the one-step input-output map; therefore the statistical bound is stated for the extended Markov operator including D. If D is fixed by the model class, the D-term below simply drops out. Let ⊂ℝdcC ^d_c be compact, let Θ⊂ℝd ^d be the compact trainable parameter set, and write Gθ(c)=(Aθ(c),Bθ(c),Cθ(c),Dθ(c)).G_θ(c)= (A_θ(c),B_θ(c),C_θ(c),D_θ(c) ). All matrix norms are operator norms unless ∥⋅∥F\|·\|_F is written. Define MB M_B :=supθ,c∥Bθ(c)∥,MC:=supθ,c∥Cθ(c)∥,MD:=supθ,c∥Dθ(c)∥, := _θ,c\|B_θ(c)\|, M_C:= _θ,c\|C_θ(c)\|, M_D:= _θ,c\|D_θ(c)\|, DΘ D_ :=supθ,θ′∈Θ‖θ−θ′‖2. := _θ,θ ∈ \|θ-θ \|_2. The tanh -bounded adapter is Lipschitz on Θ× ×C. Let LA:=supθ≠θ′,c‖Aθ(c)−Aθ′(c)‖θ−θ′‖2,LB,LC,LDbe defined analogously.L_A:= _θ≠θ ,c \|A_θ(c)-A_θ (c)\|\|θ-θ \|_2, L_B,\ L_C,\ L_D defined analogously. These constants are finite because Θ and C are compact, the adapter is a finite composition of affine maps and tanh , and ∥tanh′∥∞≤1\| \|_∞≤ 1. Let nx,nu,nyn_x,n_u,n_y denote the state, input, and output dimensions. AQS constants. Affine quadratic stability gives P=P⊤≻0P=P 0 and ε>0 >0 such that, for every θ∈Θθ∈ and c∈c , Aθ(c)⊤PAθ(c)−P⪯−εI.A_θ(c) PA_θ(c)-P - I. (24) Let λ−:=λmin(P),λ+:=λmax(P),ρ:=1−ε/λ+∈(0,1),MA:=λ+/λ−. _-:= _ (P), _+:= _ (P), ρ:= 1- / _+∈(0,1), M_A:= _+/ _-. The exogenous inputs ut∈ℝnuu_t ^n_u and output noise ξt∈ℝny _t ^n_y are centered and sub-Gaussian: ‖ut‖ψ2≤Ku,‖ξt‖ψ2≤Kξ.\|u_t\|_ _2≤ K_u, \| _t\|_ _2≤ K_ξ. For the bounded squared-loss theorem used below, assume the synthetic identification protocol has almost-sure radii ‖ut‖2≤Ru\|u_t\|_2≤ R_u and ‖ξt‖2≤Rξ\| _t\|_2≤ R_ξ. For unbounded sub-Gaussian inputs, the same proof applies after the standard truncation argument with Ru,RξR_u,R_ξ replaced by the corresponding high-probability sub-Gaussian radii. For a frozen context c, define Hθ,k(c):=Cθ(c)Aθk(c)Bθ(c),k≥0,H_θ,k(c):=C_θ(c)A_θ^k(c)B_θ(c), k≥ 0, and the extended Markov operator ℋ~θ(c):=(Dθ(c),Hθ,0(c),Hθ,1(c),…). H_θ(c):= (D_θ(c),H_θ,0(c),H_θ,1(c),… ). The well-specified observation model is yt=Dθ⋆(ct)ut+∑k=0∞Hθ⋆,k(ct)ut−k−1+ξt,θ⋆∈Θ.y_t=D_ _ (c_t)u_t+ _k=0^∞H_ _ ,k(c_t)u_t-k-1+ _t, _ ∈ . (25) The context is frozen for the local Markov operator in (25). If contexts vary inside the recurrence, the identified object is instead the time-varying product C(ct)A(ct)⋯A(ct−k+1)B(ct−k)C(c_t)A(c_t)·s A(c_t-k+1)B(c_t-k). Persistent excitation is required at the same operator level: there exists κpe>0 _ pe>0 such that, for every square-summable measurable sequence Gk(c)k≥−1\G_k(c)\_k≥-1, ‖G−1(ct)ut+∑k=0∞Gk(ct)ut−k−1‖22≥κpec[‖G−1(c)‖F2+∑k=0∞‖Gk(c)‖F2].E \|G_-1(c_t)u_t+ _k=0^∞G_k(c_t)u_t-k-1 \|_2^2≥ _ pe\,E_c [\|G_-1(c)\|_F^2+ _k=0^∞\|G_k(c)\|_F^2 ]. (26) Lemma B.1 (Coordinate-invariant equivalence class). Let V be the context index set. For two LPV realizations (A,B,C,D)(A,B,C,D) and (A′,B′,C′,D′)(A ,B ,C ,D ), define (A,B,C,D)∼(A′,B′,C′,D′)(A,B,C,D) (A ,B ,C ,D ) if and only if D(v)=D′(v)D(v)=D (v) and C(v)Ak(v)B(v)=C′(v)A′k(v)B′(v),∀v∈,k∈ℕ0.C(v)A^k(v)B(v)=C (v)A \,k(v)B (v), ∀ v ,\ k _0. Then ∼ is an equivalence relation, and each equivalence class is uniquely represented by the extended Markov operator ℋ~(v)=(D(v),C(v)Ak(v)B(v)k≥0) H(v)= (D(v),\C(v)A^k(v)B(v)\_k≥ 0 ). Proof. Reflexivity and symmetry are immediate. Transitivity follows because equality of matrices is transitive for every v and k. Hence the quotient identifies exactly those state-coordinate realizations that induce the same input-output Markov parameters. Conversely, equality of input-output maps for all finitely supported inputs identifies D(v)D(v) from the instantaneous response and identifies C(v)Ak(v)B(v)C(v)A^k(v)B(v) by applying an impulse at lag k+1k+1. Under the stochastic experiment, the same implication holds in L2L_2 by persistent excitation (26). Thus the identifiable object is the structural equivalence class represented by ℋ~(v) H(v), not a particular choice of state coordinates. ∎ Lemma B.2 (AQS implies geometric Markov decay). For every θ∈Θθ∈ , c∈c , and k≥0k≥ 0, ‖Aθk(c)‖≤MAρk,‖Hθ,k(c)‖≤MCMAMBρk.\|A_θ^k(c)\|≤ M_Aρ^k, \|H_θ,k(c)\|≤ M_CM_AM_Bρ^k. (27) The same bound holds for time-varying products Φθ,t,k=Aθ(ct)Aθ(ct−1)⋯Aθ(ct−k+1) _θ,t,k=A_θ(c_t)A_θ(c_t-1)·s A_θ(c_t-k+1). Proof. Let ‖x‖P2=x⊤Px\|x\|_P^2=x Px. From (24), ‖Aθ(c)x‖P2≤‖x‖P2−ε‖x‖22≤(1−ε/λ+)‖x‖P2=ρ2‖x‖P2.\|A_θ(c)x\|_P^2≤\|x\|_P^2- \|x\|_2^2≤ (1- / _+ )\|x\|_P^2=ρ^2\|x\|_P^2. Iterating gives ‖Aθk(c)x‖P≤ρk‖x‖P\|A_θ^\,k(c)x\|_P≤ρ^k\|x\|_P, and converting between the Euclidean and P-norms yields ‖Aθk(c)‖≤λ+/λ−ρk=MAρk\|A_θ^k(c)\|≤ _+/ _-ρ^k=M_Aρ^k. The common Lyapunov matrix P gives the identical estimate for arbitrary time-varying products. Multiplying by the uniform B- and C-bounds gives (27). ∎ Lemma B.3 (Geometric decay implies exponential β-mixing). Let Zt=(ct,ut,ξt,yt)Z_t=(c_t,u_t, _t,y_t) be the stationary process generated by (25). Suppose the exogenous process (ct,ut,ξt)(c_t,u_t, _t) is a causal process of innovations with density bounded by Kp<∞K_p<∞ and absolute-regularity coefficients βexo(q)≤b0e−a0q _ exo(q)≤ b_0e^-a_0q; the independent-exogenous case is obtained by taking b0=0b_0=0. Then ZtZ_t is exponentially β-mixing: there exist constants bβ<∞b_β<∞ and ρβ∈(0,1) _β∈(0,1), depending only on (MA,MB,MC,MD,ρ,Ru,Rξ,b0,a0,Kp,nu,ny,dc)(M_A,M_B,M_C,M_D,ρ,R_u,R_ξ,b_0,a_0,K_p,n_u,n_y,d_c), such that βZ(q)≤bβρβq,q≥1. _Z(q)≤ b_β _β^q, q≥ 1. (28) One admissible choice is ρβ:=maxρ1/2,e−a0/2,bβ:=b0+CPT(Kp,nu,ny,dc)(1+MCMAMBRu1−ρ+MDRu+Rξ) _β:= \ρ^1/2,e^-a_0/2\, b_β:=b_0+C_ PT(K_p,n_u,n_y,d_c) (1+ M_CM_AM_BR_u1-ρ+M_DR_u+R_ξ ) where CPTC_ PT is the finite constant in the Pham–Tran absolute-regularity theorem for Volterra/Bernoulli-shift processes [29]. Proof. Let yt(q)y_t^(q) be the output obtained by truncating (25) to inputs no older than q lags. By Lemma B.2, ‖yt−yt(q)‖2≤∑k>q‖Hθ⋆,k(ct)‖‖ut−k−1‖2≤MCMAMBRu1−ρρq+1.\|y_t-y_t^(q)\|_2≤ _k>q\|H_ _ ,k(c_t)\|\,\|u_t-k-1\|_2≤ M_CM_AM_BR_u1-ρ^q+1. (29) Thus ZtZ_t is a causal Bernoulli-shift, equivalently a Volterra process, whose remote-past influence is exponentially summable. The theorem of Pham and Tran [29] states that a Volterra process with bounded-density innovations and summable kernel tails is absolutely regular, with mixing coefficients controlled by the innovation mixing coefficient plus the kernel tail. Combining βexo(q)≤b0e−a0q _ exo(q)≤ b_0e^-a_0q with (29) and enlarging constants gives (28). ∎ Lemma B.4 (Dependent Rademacher complexity). For zt=(ct,ut−jj≥0,yt)z_t=(c_t,\u_t-j\_j≥ 0,y_t), define fθ(zt)=Dθ(ct)ut+∑k=0∞Hθ,k(ct)ut−k−1,ℓθ(zt):=‖yt−fθ(zt)‖22.f_θ(z_t)=D_θ(c_t)u_t+ _k=0^∞H_θ,k(c_t)u_t-k-1, _θ(z_t):=\|y_t-f_θ(z_t)\|_2^2. Let ℒ:=ℓθ:θ∈ΘL:=\ _θ:θ∈ \. Define Bf:=(MD+MCMAMB1−ρ)Ru,By:=Bf+Rξ,Bℓ:=4By2,B_f:= (M_D+ M_CM_AM_B1-ρ )R_u, B_y:=B_f+R_ξ, B_ :=4B_y^2, and Lf:=Ru[LD+LCMAMB+MCMALB1−ρ+MCMBLAMA2(1−ρ)2].L_f:=R_u [L_D+ L_CM_AM_B+M_CM_AL_B1-ρ+ M_CM_BL_AM_A^2(1-ρ)^2 ]. Finally let Γβ _β :=1+2∑q=1∞βZ(q)≤1+2bβρβ1−ρβ, :=1+2 _q=1^∞ _Z(q)≤ 1+ 2b_β _β1- _β, Jd J_d :=∫01dlog(1+3LfDΘBfs)ds. := _0^1 d (1+ 3L_fD_ B_fs )\,ds. Then the empirical Rademacher complexity of the squared-loss class satisfies ℜNβ(ℒ)≤CradN, R_N^β(L)≤ C_ rad N, (30) where Crad:=96ByBfnyΓβJd.C_ rad:=96\,B_yB_f n_y _β\,J_d. Proof. The geometric series in Lemma B.2 gives ‖fθ(zt)‖2≤Bf\|f_θ(z_t)\|_2≤ B_f and ‖yt‖2≤By\|y_t\|_2≤ B_y, hence 0≤ℓθ≤Bℓ0≤ _θ≤ B_ . For any θ,θ′θ,θ , ‖Hθ,k(c)−Hθ′,k(c)‖ \|H_θ,k(c)-H_θ ,k(c)\| ≤(LCMAMBρk+MCMALBρkCLOSE ≤ (L_CM_AM_Bρ^k+M_CM_AL_Bρ^k OPEN+MCMBLAMA2kρk−1)‖θ−θ′‖2, +M_CM_BL_AM_A^2kρ^k-1 )\|θ-θ \|_2, and summing over k≥0k≥ 0, adding the D-term, and multiplying by RuR_u gives ‖fθ(zt)−fθ′(zt)‖2≤Lf‖θ−θ′‖2.\|f_θ(z_t)-f_θ (z_t)\|_2≤ L_f\|θ-θ \|_2. Thus the covering number of each output coordinate satisfies log(ϵ,fθ,j,L2)≤dlog(1+3LfDΘϵ). (ε,\f_θ,j\,L_2)≤ d (1+ 3L_fD_ ε ). Dudley’s entropy integral gives ℜN(fθ,j)≤12BfNJd. R_N(\f_θ,j\)≤ 12B_f NJ_d. Summing over nyn_y output coordinates and using Cauchy–Schwarz gives the ny n_y factor. Since ℓθ(z)=‖y−fθ(z)‖22 _θ(z)=\|y-f_θ(z)\|_2^2 is 4By4B_y-Lipschitz in fθ(z)f_θ(z) on the bounded range, Talagrand’s contraction principle [22, 23] transfers the predictor bound to the squared-loss class. For the dependent sequence, the independent-block coupling of Yu and the Mohri–Rostamizadeh blocking lemma for β-mixing sequences inflate the independent complexity by at most Γβ _β [30, 31]. This yields (30). ∎ Proposition B.1 (Excess-risk and operator convergence). Let θ θ be an empirical risk minimizer over Θ : θ^∈argminθ∈Θℛ^N(θ),ℛ^N(θ):=1N∑t=1Nℓθ(zt), θ∈ _θ∈ R_N(θ), R_N(θ):= 1N _t=1^N _θ(z_t), and let ℛ(θ):=ℓθ(Zt)R(θ):=E _θ(Z_t). For any δ∈(0,1)δ∈(0,1), with probability at least 1−δ1-δ, ℛ(θ^)−ℛ(θ⋆)≤Cex(δ)N,R( θ)-R( _ )≤ C_ ex(δ) N, (31) where Cex(δ):=4Crad+2Bℓ2Γβlog(2/δ).C_ ex(δ):=4C_ rad+2B_ 2 _β (2/δ). Moreover, the squared extended Markov-operator error obeys c[‖Dθ^(c)−Dθ⋆(c)‖F2+∑k=0∞‖Hθ^,k(c)−Hθ⋆,k(c)‖F2]≤Cex(δ)κpeN.E_c [\|D_ θ(c)-D_ _ (c)\|_F^2+ _k=0^∞\|H_ θ,k(c)-H_ _ ,k(c)\|_F^2 ]≤ C_ ex(δ) _ pe N. (32) Proof. For bounded losses generated by an exponentially β-mixing process, the uniform convergence theorem for absolutely regular sequences gives [31, 32] supθ∈Θ|ℛ(θ)−ℛ^N(θ)|≤2ℜNβ(ℒ)+Bℓ2Γβlog(2/δ)N _θ∈ |R(θ)- R_N(θ) |≤ 2 R_N^β(L)+B_ 2 _β (2/δ)N with probability at least 1−δ1-δ. Substituting Lemma B.4 gives the right-hand side ΔN(δ)≤2Crad+Bℓ2Γβlog(2/δ)N. _N(δ)≤ 2C_ rad+B_ 2 _β (2/δ) N. Since θ θ minimizes the empirical risk and θ⋆∈Θ _ ∈ , ℛ(θ^)−ℛ(θ⋆)≤2ΔN(δ)≤Cex(δ)N,R( θ)-R( _ )≤ 2 _N(δ)≤ C_ ex(δ) N, which proves (31). It remains to connect the statistical risk to the coordinate-invariant operator. Because the noise is centered and independent of the inputs, and the model is well specified, the noise variance cancels in the excess risk: ℛ(θ)−ℛ(θ⋆)=‖ΔDθ(ct)ut+∑k=0∞ΔHθ,k(ct)ut−k−1‖22.R(θ)-R( _ )=E \| D_θ(c_t)u_t+ _k=0^∞ H_θ,k(c_t)u_t-k-1 \|_2^2. Here ΔDθ=Dθ−Dθ⋆ D_θ=D_θ-D_ _ and ΔHθ,k=Hθ,k−Hθ⋆,k H_θ,k=H_θ,k-H_ _ ,k. Applying persistent excitation (26) with G−1=ΔDθG_-1= D_θ and Gk=ΔHθ,kG_k= H_θ,k yields ℛ(θ)−ℛ(θ⋆)≥κpec[‖ΔDθ(c)‖F2+∑k=0∞‖ΔHθ,k(c)‖F2].R(θ)-R( _ )≥ _ pe\,E_c [\| D_θ(c)\|_F^2+ _k=0^∞\| H_θ,k(c)\|_F^2 ]. Taking θ=θ^θ= θ and using (31) proves (32). ∎ Combining Lemma B.1 with Proposition B.1 proves Proposition 4.1. The internal matrices are identifiable only modulo the structural equivalence relation ∼ , while the observable Markov operator ℋH is identifiable. Under squared loss, both the excess prediction risk and the squared extended Markov-operator error are bounded by C/NC/ N with high probability, where C is an explicit function of the AQS certificate (P,ε)(P, ), representation bounds (MA,MB,MC,MD)(M_A,M_B,M_C,M_D), excitation and tail constants (Ku,Kξ,Ru,Rξ,κpe)(K_u,K_ξ,R_u,R_ξ, _ pe), the parameter dimension d, the tanh -parameterization Lipschitz constants (LA,LB,LC,LD)(L_A,L_B,L_C,L_D), the mixing constants (bβ,ρβ)(b_β, _β), and the confidence level δ, but not of N. The unsquared Frobenius norm is controlled by taking a square root of (32); the rate stated in Proposition 4.1 is the squared-loss estimation rate. Appendix C Experiment Details C.1 Training Details Training adapts the frozen Hydra-BERT backbone to the masked diffusion objective on an internal Hydra corpus constructed from the source mixture in Table 5. The reported train and validation portions contain 4,295,110,656 tokens in total: 3,537,149,952 training tokens and 757,960,704 validation tokens. The validation split is used for monitoring and reporting the held-out losses in the experiments, while all model updates use only the training split. Table 5: Source mixture used to construct the Hydra train/validation corpus. Percentages are the relative shard-pool weights of the reported training runs and normalize at sampling. Source Ratio (%) common-pile StackExchange 10.00 common-pile LibreTexts 1.67 common-pile YouTube 22.50 common-pile PubMed 7.50 common-pile C 10.00 common-pile Project Gutenberg 5.00 common-pile arXiv papers 7.50 common-pile News 1.67 common-pile DOAB 15.00 common-pile Pressbooks 1.67 WikiText-103 raw 2.50 common-pile Data Provenance Initiative 8.33 common-pile Wikimedia 8.33 common-pile arXiv abstracts 3.33 For all reported experiments, only the MaRK adapter and timestep-embedding parameter groups are trainable. The Hydra backbone parameters and masked-LM classification head remain fixed, isolating the effect of dynamic operator modulation. The objective is the masked-token diffusion negative log-likelihood with CART weighting enabled using parameter p=0.45p=0.45. Training uses Distributed Data Parallel (DDP) on two NVIDIA GH200 GPUs with bfloat16 precision, gradient checkpointing, high-precision matrix multiplication, and gradient clipping at norm 1.01.0. Table 6: Training settings shared by all three MaRK kernels. Values are taken from the training configurations and the corresponding optimization implementation. Setting Value Trainable parameters MaRK adapters and timestep embeddings only Backbone Frozen 23-layer Hydra-BERT SSM Sequence length 4096 tokens Epochs 3 Optimizer Fused AdamW, β=(0.9,0.999)β=(0.9,0.999), weight decay 0.010.01 MaRK learning rate 6.0×10−46.0× 10^-4 Learning-rate schedule Step-wise linear warmup, then cosine decay Minimum LR ratio 0.1 Precision bfloat16 Training strategy Distributed Data Parallel Workers 8 Matrix multiply precision high Hydra hidden size / vocabulary 768 / 30,522 SSM state / convolution / head dim 64 / 7 / 64 Expansion / chunk size 2 / 256 The variants are initialized from their corresponding kernel-specific checkpoints. Apart from the kernel geometry, all reported variants use the same batch size, warmup length, total training horizon, optimizer, and base learning rate. The 6.0×10−46.0× 10^-4 learning rate is linearly warmed up for 720 steps and then cosine-decayed over the remaining steps of the 14,400-step schedule to a minimum value of 0.1×0.1× the base rate. The kernel-specific entries in Table 7 specify only the adapter geometry: Chebyshev uses a degree-5 polynomial basis, while DCT uses a 256-point basis grid with 8 retained frequencies. Table 7: Kernel-specific training settings. Shared kernel dimensions are rank r=2r=2, adapter MLP dimension 256, and timestep embedding dimension 128. Kernel Batch size Warmup / total steps Kernel-specific parameters Hypernet 112 720 / 14,400 – Chebyshev 112 720 / 14,400 Polynomial degree =5=5 DCT 112 720 / 14,400 Timepoints =256=256, frequencies =8=8 C.2 Modulation Bounds and CART Settings The bounded shifts and scales of Section 2.1 use the magnitudes in Table 8. Diagonal parameters use βΘtanh(⋅) _ (·), with the A-shift bound capped at βA≤0.5 _A≤ 0.5 and the remaining magnitudes learned; dense parameters use αΦ _ and βΦ _ on the low-rank factors. Every magnitude is initialized to 0.20.2. The kernel dimensions (low-rank size, adapter MLP, timestep embedding, Chebyshev degree, and DCT frequency budget) are listed in Table 7. Table 8: Modulation bounds and CART settings, identical across the three kernels. Parameter Value A-shift bound βA _A ≤0.5≤ 0.5 Scale and shift initial magnitude α,βα,β 0.20.2 (learned) CART p / scale 0.45 / 1.0 The A-shift bound is the only constraint that enters the stability certificate, so it is capped while the remaining magnitudes are learned freely; B, C, and D never enter the Lyapunov condition because they are not raised to powers. Because A<0A<0 and Δ>0 >0, every discrete eigenvalue stays in (0,1)(0,1) under any finite bound, so the modulated recurrence is Schur-stable for each individual parameter value and the Markov sequence always decays. The finite bound buys a uniform positive margin rather than stability itself: the worst case sits at the vertex δ=−βAδ=- _A, and as βA _A grows, that eigenvalue approaches 11 from below while the margin ε=1−maxiA¯i2 =1- _i A_i^2 shrinks toward zero. An unbounded shift would therefore not destabilize any single realization, but it would remove the uniform guarantee that Appendix B relies on, which is why the bound is kept finite. The bounds themselves were set conservatively to keep the modulated operators near the frozen backbone, and the CART sharpness p=0.45p=0.45 follows Dream. C.3 Evaluation Protocol and Statistics Every language-model number in this paper is a CART-weighted validation loss, read as weighted_nll from the evaluation harness; perplexity is only the exponential of that quantity and is not the reported metric. Each cell is a mean over 1010 seeds for the accuracy tables (Tables 1 and 2) and over 1515 seeds for the inference table (Table 3), with seed i fixed at 42+100i42+100i. Intervals are standard-error-of-the-mean bands over seeds (sample standard deviation with ddof=1ddof=1 divided by n n): Tables 1 and 2 report the normal approximation (mean ± 1.96×±\,1.96×SEM over 1010 seeds), while the inference table (Table 3) reports the Student-t interval with 1414 degrees of freedom over 1515 seeds. The synthetic LPV curves in Figure 2 use the same seed and interval convention. Validation consumes the global random-number stream on every batch, sampling both the diffusion timestep and the masking pattern, and building the model also draws from it because Hydra initializes its discretization bias randomly. The harness therefore captures the RNG state immediately after the model is built and restores it before each dataset, so the CART suite reproduces the leave-one-out draws and the WikiText loss matches the leave-one-out full row to the bfloat16 run-to-run noise. Because the largest validation split (arXiv) contains tens of thousands of packed sequences at length 4096, its evaluation is capped at 500 validation batches, while smaller datasets are evaluated in full. C.4 AQS Certificate Evaluation We evaluate the AQS certificate by directly checking discrete-time stability margins implied by the ZOH discretization for the learned diagonal recurrence. Using the identity Lyapunov function P=IP=I, this yields a closed-form certificate that avoids SDP solving. We load the pretrained diagonal log-parameters Alog(ℓ)A_ ^( ) from the 23-layer Hydra checkpoint and evaluate the worst-case (least stable) parameter vertex by monotonicity of the map Alog↦A¯=exp(−exp(Alog)Δ)A_ A= (- (A_ ) ) under the bounded shift Alog±βA_ ± _A, using Δmin=0.001 _ =0.001: Algorithm 4 Constructive AQS Certificate via P=IP=I Lyapunov Function 1: Input: Pretrained Alog(ℓ)A_ ^( ) for layers ℓ=1,…,23 =1,…,23; bound βA=0.5 _A=0.5; Δmin=0.001 _ =0.001 2: Set P←IP← I ⊳ Identity matrix (valid Lyapunov function for diagonal A with ρ<1ρ<1) 3: for each layer ℓ , each head i do 4: Alog,iworst←Alog,i(ℓ)−βA_ ,i^worst← A_ ,i^( )- _A ⊳ Worst-case vertex (least decay) 5: Ai←−exp(Alog,iworst)A_i←- (A_ ,i^worst) ⊳ Strictly negative 6: a¯i←exp(Ai⋅Δmin) a_i← (A_i· _ ) ⊳ Discrete eigenvalue ∈(0,1)∈(0,1) 7: end for 8: ρmax(ℓ)←maxi|a¯i| _ ^( )← _i| a_i| ⊳ Spectral radius at worst-case vertex 9: εℓ←1−(ρmax(ℓ))2 _ ← 1-( _ ^( ))^2 ⊳ Lyapunov stability margin 10: Verify: εℓ>0 _ >0 for all ℓ ⊳ A¯⊤IA¯−I≺0 A I\, A-I 0 (neg. definite) 11: Return P=IP=I, εℓ=123\ _ \_ =1^23 ⊳ AQS Certificate with per-layer margins The procedure is conservative in three ways: it fixes the Lyapunov function to P=IP=I, evaluates the least stable vertex δ=−βAδ=- _A, and uses the smallest discretization Δmin=0.001 _ =0.001, which places the discrete eigenvalue closest to one. Because βA _A is shared across the bounded kernels, the certificate depends only on the pretrained AlogA_ and not on the adapter geometry, so a single pass over the backbone covers all three variants. The reported margins are the distance of the worst-case eigenvalue from the unit circle. C.5 Synthetic LPV Identifiability Algorithms For identifiability tests, both ground-truth generator models and target estimators use diagonal LPV-SSMs with state dimension dstate=4d_state=4, factorization rank r=2r=2, sequence length L=20L=20, and datasets N∈50,100,200,400,800,1600N∈\50,100,200,400,800,1600\ over 15 random seeds. Inputs are sampled as u∼(0,1)u (0,1) and diffusion coordinates as τ∼(0,1)τ (0,1); scalar outputs are generated by the ground-truth system and corrupted with Gaussian noise σ=0.1σ=0.1. Chebyshev and DCT use a fixed readout C=C=1 to remove final rotation/scaling ambiguity, while the Hypernet synthetic kernel also modulates C(τ)C(τ) with bounded low-rank scale/shift terms. This choice is without loss of generality for scalar readouts: any nonzero scalar C can be absorbed by a linear reparameterization of the latent state, so fixing C=C=1 isolates operator recovery rather than an arbitrary output scaling. The continuous matrices returned by Algorithms 5–7 are rolled out with the same scalar ZOH discretization used in the scripts: A¯(τ)=exp(A(τ)Δt),B¯(τ)=A¯(τ)−1A(τ)+10−7⊙B(τ),Δt=0.1. A(τ)= (A(τ) t), B(τ)= A(τ)-1A(τ)+10^-7 B(τ), t=0.1. The synthetic sequence is then generated by xs=A¯(τs)⊙xs−1+B¯(τs)us,ys=C(τs)⊤xs.x_s= A( _s) x_s-1+ B( _s)u_s, y_s=C( _s) x_s. Estimators are initialized from noisy copies of the ground-truth parameters and trained with Mean-Squared Error (MSE). Recovery is measured by frozen-time Markov operator error on a 200-point τ grid over lags k=0,…,14k=0,…,14, using Mk(τ)=C(τ)⊤(Ak¯(τ)⊙B¯(τ))M_k(τ)=C(τ) ( A^k(τ) B(τ)); for Chebyshev and DCT, C(τ)=C(τ)=1 is constant. The Hypernet generator is the most flexible synthetic parameterization: a sinusoidal embedding of τ is mapped by a four-layer SiLU MLP to bounded shifts of the recurrence, input map, and readout. The negative exponential form enforces A(τ)<0A(τ)<0 componentwise, and sigmoid/tanh bounds keep the modulation amplitude controlled. Unlike the basis-kernel variants below, Hypernet returns a time-varying C(τ)C(τ) in addition to A(τ)A(τ) and B(τ)B(τ). Algorithm 5 Hypernet Synthetic LPV-SSM 1: Input: Diffusion coordinate τ∈[0,1]τ∈[0,1] 2: Parameters: Alog_baseA_ \_base, low-rank base Ubase,VbaseU_base,V_base, scalars βAparam,αB,βB,αC,βC _A^param, _B, _B, _C, _C, four-layer SiLU MLP, linear heads 3: cond(τ)←[sin(2πmτ),cos(2πmτ)]m=063cond(τ)←[ (2π mτ), (2π mτ)]_m=0^63 ⊳ 128-dimensional conditioning signal 4: h←MLP4×SiLU(cond(τ))h _4×SiLU(cond(τ)) 5: uvflat(q)←vec(UqVq)1:dstateuvflat(q) (U_qV_q)_1:d_state, where q is split into rank-r factors (Uq,Vq)(U_q,V_q) ⊳ matches the low-rank head flattening 6: βAeff←0.5σ(βAparam) _A^eff← 0.5\,σ( _A^param), ashift←WAha_shift← W_Ah ⊳ bounds the recurrence modulation 7: A(τ)←−exp(Alog_base+βAefftanh(ashift))A(τ)←- (A_ \_base+ _A^eff (a_shift)) 8: bbase←vec(UbaseVbase)1:dstateb_base (U_baseV_base)_1:d_state 9: bscale(τ)←1+αBtanh(uvflat(WB,scaleh))b_scale(τ)← 1+ _B (uvflat(W_B,scaleh)), bshift(τ)←βBtanh(uvflat(WB,shifth))b_shift(τ)← _B (uvflat(W_B,shifth)) ⊳ bounded scale/shift for B 10: B(τ)←bbase⊙bscale(τ)+bshift(τ)B(τ)← b_base b_scale(τ)+b_shift(τ) 11: cscale(τ)←1+αCtanh(uvflat(WC,scaleh))c_scale(τ)← 1+ _C (uvflat(W_C,scaleh)), cshift(τ)←βCtanh(uvflat(WC,shifth))c_shift(τ)← _C (uvflat(W_C,shifth)) 12: C(τ)←⊙cscale(τ)+cshift(τ)C(τ) 1 c_scale(τ)+c_shift(τ) ⊳ learned readout for the Hypernet variant 13: Return A(τ),B(τ),C(τ)A(τ),B(τ),C(τ) Hypernet is the only synthetic generator with a learned, time-varying readout, so its recovery error also reflects the readout path and not just the recurrence and input maps. It is therefore the hardest of the three synthetic targets, and the basis kernels are compared against it rather than against an easier diagonal reference. The Chebyshev generator uses a fixed polynomial basis after mapping τ to z∈[−1,1]z∈[-1,1]. All time variation enters through bounded Chebyshev shifts of the log-recurrence and of the low-rank factors defining B(τ)B(τ). The readout is fixed to C=C=1, so the reported recovery error measures the induced Markov operator rather than an arbitrary rescaling of the output map. Algorithm 6 Chebyshev Synthetic LPV-SSM 1: Input: Diffusion coordinate τ∈[0,1]τ∈[0,1] 2: Parameters: Alog_baseA_ \_base, UbaseU_base, VbaseV_base, coefficients Acoeffs,Ucoeffs,VcoeffsA_coeffs,U_coeffs,V_coeffs, bounds βA,βB _A, _B, basis count Kbasis=3K_basis=3 3: z←2τ−1z← 2τ-1 ⊳ Chebyshev domain 4: T0(z)←1T_0(z)← 1, T1(z)←zT_1(z)← z 5: for k=2k=2 to Kbasis−1K_basis-1 do 6: Tk(z)←2z⋅Tk−1(z)−Tk−2(z)T_k(z)← 2z· T_k-1(z)-T_k-2(z) 7: end for 8: Ashift←∑k=0Kbasis−1Tk(z)Acoeffs(k)A_shift← _k=0^K_basis-1T_k(z)A_coeffs^(k) ⊳ basis modulation of the log recurrence 9: Ushift←∑k=0Kbasis−1Tk(z)Ucoeffs(k)U_shift← _k=0^K_basis-1T_k(z)U_coeffs^(k), Vshift←∑k=0Kbasis−1Tk(z)Vcoeffs(k)V_shift← _k=0^K_basis-1T_k(z)V_coeffs^(k) 10: A(τ)←−exp(Alog_base+βAtanh(Ashift))A(τ)←- (A_ \_base+ _A (A_shift)) ⊳ negative diagonal continuous dynamics 11: U(τ)←Ubase+βBtanh(Ushift)U(τ)← U_base+ _B (U_shift), V(τ)←Vbase+βBtanh(Vshift)V(τ)← V_base+ _B (V_shift) 12: B(τ)←U(τ)V(τ)B(τ)← U(τ)V(τ), C←C 1 ⊳ fixed readout removes scaling ambiguity 13: Return A(τ),B(τ),CA(τ),B(τ),C The synthetic Chebyshev generator uses a degree-3 basis to keep the target low-dimensional, whereas the language-model Chebyshev adapter uses degree 5 (Table 7). The synthetic target is therefore a smoother, lower-capacity member of the same polynomial family. The DCT generator replaces the Chebyshev basis with cosine modes and a learned monotone spectral decay. The normalization constants match the Ltime=256L_time=256 DCT-I basis used in the experiment, while tanh-bounded shifts keep the recurrence and input map controlled. As in the Chebyshev setting, C=C=1 is fixed so the comparison focuses on the recovered operator. Algorithm 7 Discrete Cosine Transform (DCT) Synthetic LPV-SSM 1: Input: Diffusion coordinate τ∈[0,1]τ∈[0,1] 2: Parameters: Alog_baseA_ \_base, UbaseU_base, VbaseV_base, spectral coefficients Aspec,Uspec,VspecA_spec,U_spec,V_spec, bounds βA,βB _A, _B, bandwidth α, nfreqs=8n_freqs=8, Ltime=256L_time=256 3: for k=0k=0 to nfreqs−1n_freqs-1 do 4: wk←(k+1)−(softplus(α)+10−4)w_k←(k+1)^-(softplus(α)+10^-4) ⊳ learned spectral attenuation 5: ck←1/Ltimec_k← 1/ L_time if k=0k=0, else 2/Ltime 2/L_time ⊳ DCT-I normalization 6: bk(τ)←ckcos(πkτ)b_k(τ)← c_k (π kτ) 7: end for 8: Ashift←∑k=0nfreqs−1bk(τ)wkAspec(k)A_shift← _k=0^n_freqs-1b_k(τ)\,w_k\,A_spec^(k) ⊳ attenuated spectral recurrence shift 9: Ushift←∑k=0nfreqs−1bk(τ)wkUspec(k)U_shift← _k=0^n_freqs-1b_k(τ)\,w_k\,U_spec^(k), Vshift←∑k=0nfreqs−1bk(τ)wkVspec(k)V_shift← _k=0^n_freqs-1b_k(τ)\,w_k\,V_spec^(k) 10: A(τ)←−exp(Alog_base+βAtanh(Ashift))A(τ)←- (A_ \_base+ _A (A_shift)) ⊳ negative diagonal continuous dynamics 11: U(τ)←Ubase+βBtanh(Ushift)U(τ)← U_base+ _B (U_shift), V(τ)←Vbase+βBtanh(Vshift)V(τ)← V_base+ _B (V_shift) 12: B(τ)←U(τ)V(τ)B(τ)← U(τ)V(τ), C←C 1 ⊳ fixed readout isolates operator recovery 13: Return A(τ),B(τ),CA(τ),B(τ),C The synthetic DCT generator keeps the basis length and frequency count of the language-model adapter so that the target uses the same bandlimited parameterization it is meant to probe. The learned attenuation only reshapes how quickly the high-frequency coefficients decay, so a miss in recovery is attributable to the bandlimited basis rather than to a mismatch in basis definition.