Modeling Dynamic Mixtures of Time-Delay Systems from Streaming Time Series

# Modeling Dynamic Mixtures of Time-Delay Systems from Streaming Time Series
Source: [https://arxiv.org/html/2605.26191](https://arxiv.org/html/2605.26191)
Ren FujiwaraYasuko MatsubaraYasushi Sakurai SANKEN, The University of Osaka, Japan \{r\-fujiwr,yasuko,yasushi\}@sanken\.osaka\-u\.ac\.jp

###### Abstract

This research addresses the problem of adaptive modeling in time\-series data streams with clear input\-output relationships\. This problem is challenging because rapid system changes \(regime shifts\) caused by environmental factors or input delay changes degrade model performance, and the trade\-off among accuracy, robustness, and memory usage arises when using multiple small models for each time\-series pattern\. To address these issues, this paper presents an online framework/method that treats streaming time series as dynamic mixtures of time\-delay systems\. This framework maintains robustness of model tracking and reduces memory usage by summarizing past regimes using a fixed\-length representation that captures both the system dynamics and input\-output delays\. Concretely, this approach constructs a summary system tensor using the system’s Markov parameter series, capturing both dynamic behavior and delay characteristics\. If necessary, a tensor decomposition algorithm extracts relevant past models from the tensor and helps select the system that best fits the current regime\. This method enables rapid adaptation to environmental changes and is computationally efficient\. Tests on real datasets show thatDelayMixconsistently outperforms other methods, achieving superior forecast accuracy and faster adaptation to delays, especially for highly non\-stationary data\.

## 1Introduction

In recent years, the importance of modeling time\-series data streams with explicit input\-output relationships and modeling based on such data streams has grown\. In particular, in real\-world systems such as ship controlli2022hybrid, industrial robotslu2021gan, data center coolingdatacentor, diesel engine controlengine, and autonomous drivingdrivectr, modeling failures can easily manifest as failures in safe and efficient operation, necessitating the development of methods for stable and efficient modeling on time\-series data streams\.

However, modeling on streams simultaneously faces the challenges of data nonstationarity and computational resource constraints, making it difficult to directly apply models and algorithms designed for offline use\. Specifically, the following two problems arise: \(1\) Changes in the underlying system often cause abrupt model changes, or regime shifts, rapidly degrading the performance of trained modelsRegimeShift;RegimeShift2;RegimeShiftNature\. Furthermore, this challenge is further complicated by the unknown input delay inherent in systems with input\-output relationshipsARIMAX;MPC;kLinReg;QLag;SINDyc;DMDc;ReLiNet\. This means that changes in the model’s input delays must be simultaneously tracked\. \(2\) Time\-series data streams have no temporal boundaries \(i\.e\., they are semi\-infinitely long\), making model estimation using all past data prone to design errors\. Furthermore, models designed using multiple small models for each time\-series pattern suffer from a trade\-off between robustness of accuracy and memory usage\. While retaining more models stabilizes accuracy, having too many models increases the model search cost and memory computational complexity without bound\.

For example, current stream prediction models that consider regime shifts offer scalable algorithms in terms of model adaptability and computational time on the streamRegimeCast;regime\_other;ModePlait;MicroAdapt\. However, these methods monotonically increase the number of models, thereby increasing not only memory usage but also the computational time required for model switching\. Furthermore, these methods do not consider exogenous variables\. In other words, these methods lack design principles for tracking changes in input delay, making them difficult to directly apply to our problem\. Recent neural approaches, such as ReLiNetReLiNetand TimeXerTimeXer, incorporate exogenous signals through switched linear systems and Transformer\-based attention, respectively, achieving competitive results\. However, these models generally do not model time delays explicitly and tend to struggle in streaming settings due to high computational cost and frequent retraining\. Therefore, this paper aims to develop a method that simultaneously satisfies both adaptability to sudden environmental changes and scalable estimation\. The research question in this paper is:

Can we develop a method that achieves high\-accuracy modeling of time\-series data streams with explicit input\-output relationships, while tracking rapid environmental changes and maintaining high computational efficiency?

We proposeDelayMix, an online framework that resolves the conflict between structural identification and streaming adaptation\. This framework maintains high model tracking performance while avoiding increased memory usage by retaining model sets used in past regimes as mixtures of fixed\-length summarized representations of the system’s dynamic behavior and input\-output delay structure\. Specifically, we focus on the fact that the system’s Markov parameter series integrates its dynamic behavior and input\-output delay structure, and we construct a system tensor that summarizes them\. Then, if necessary, we use a tensor decomposition algorithm to extract useful past models from the system tensor and adaptively select the system that best fits the current regime\. This approach achieves both tracking rapid environmental changes and high computational efficiency\.

In summary,DelayMixis an online algorithm for modeling streaming time series as mixtures of time\-delay systems\. Our main contributions are as follows:

- •Streaming Modeling of Time\-Delay Systems:We propose a new formulation that leverages Markov parameter sequences to jointly learn unknown delays and system dynamics in streaming environments\.
- •Scalability with Resource Constraints:The computation time ofDelayMixis independent of the time series length, thanks to incremental model updating, and thus is a faster algorithm than its competitors\.
- •Robustness Against Structural Drift:Experiments on real\-world datasets demonstrate that our method significantly outperforms state\-of\-the\-art baselines in predictive accuracy, especially in environments with strong nonstationarity\.

## 2Proposed Method

In this paper, we study time series arising from mixtures of time\-delay MIMO systems, where the delays are unknown and may differ across systems\. Our goal is to recover the dynamics of each system from one observed sequence and use them for online multi\-step forecasting\. Rather than estimating delays directly, we estimate Markov parameters for each regime to summarize input–output behavior\. For each regime, we then create a delay\-free input–output equivalent state\-space model, which we use for inference and prediction\.

This section introduces the proposed model and the algorithm\. First, in section[2\.1](https://arxiv.org/html/2605.26191#S2.SS1), we define our model\. Then we formalize the problem definition in Problem[1](https://arxiv.org/html/2605.26191#Thmproblem1)\. We then present the proposed online algorithmDelayMix, which includes sub\-algorithms named DynamicMomentCollection, ModelDatabaseAdaption, and FuturePrediction\. The theoretical justification for our method is given in Section[3](https://arxiv.org/html/2605.26191#S3)\.

### 2\.1Mixtures of Time\-Delay Systems

In practice, time\-series data rarely come from just one stationary dynamical system\. Instead, the underlying dynamics often switch between several regimes, depending on hidden conditions such as the operating mode, environment, or user behavior\. To capture this non\-stationarity, we model the data as coming from a mixture of time\-delay systems\.

LetRRdenote the number of regimes, and letrt∈\{1,…,R\}r\_\{t\}\\in\\\{1,\\dots,R\\\}denote the \(unobserved\) regime index at timett\. For each regimei∈\{1,…,R\}i\\in\\\{1,\\dots,R\\\}, we introduce a latent state𝐳i​\(t\)∈ℝk\\mathbf\{z\}\_\{i\}\(t\)\\in\\mathbb\{R\}^\{k\}and system matrices𝐀i∈ℝk×k,𝐁i∈ℝk×dc,𝐂i∈ℝd×k\.\\mathbf\{A\}\_\{i\}\\in\\mathbb\{R\}^\{k\\times k\},\\mathbf\{B\}\_\{i\}\\in\\mathbb\{R\}^\{k\\times d\_\{c\}\},\\mathbf\{C\}\_\{i\}\\in\\mathbb\{R\}^\{d\\times k\}\.We also include an integer\-valued delayτi≥0\\tau\_\{i\}\\geq 0in the state update\. Conditioned on being in regimert=ir\_\{t\}=iat timett, the dynamics evolve as

𝐳i​\(t\+1\)\\displaystyle\\mathbf\{z\}\_\{i\}\(t\+1\)=𝐀i​𝐳i​\(t\)\+𝐁i​𝐮​\(t−τi\),\\displaystyle=\\mathbf\{A\}\_\{i\}\\mathbf\{z\}\_\{i\}\(t\)\+\\mathbf\{B\}\_\{i\}\\mathbf\{u\}\(t\-\\tau\_\{i\}\),\(1\)𝐱​\(t\)\\displaystyle\\mathbf\{x\}\(t\)=𝐂i​𝐳i​\(t\),\\displaystyle=\\mathbf\{C\}\_\{i\}\\mathbf\{z\}\_\{i\}\(t\),\(2\)where𝐮​\(t\)∈ℝdc\\mathbf\{u\}\(t\)\\in\\mathbb\{R\}^\{d\_\{c\}\}is the observed exogenous input and𝐱​\(t\)∈ℝd\\mathbf\{x\}\(t\)\\in\\mathbb\{R\}^\{d\}is the observed output\. The input𝐮​\(t\)\\mathbf\{u\}\(t\)is observed at every time step, but the regime indexrtr\_\{t\}and the latent states𝐳i​\(t\)\\mathbf\{z\}\_\{i\}\(t\)are not directly observed\. We use𝐱^​\(t\)\\mathbf\{\\hat\{x\}\}\(t\)for model predictions \(for example, predictions produced by filtering\) when needed later\.

Regime\-specific Markov parameters and mixture view\.For each regime, the time\-delay system \([1](https://arxiv.org/html/2605.26191#S2.E1)\)–\([2](https://arxiv.org/html/2605.26191#S2.E2)\) yields the following equivalent Markov parameters \(impulse response\):

\{𝐠1\(i\),𝐠2\(i\),…\},𝐠ℓ\(i\)=𝐂i​\(𝐀i\)ℓ−1​𝐁i\.\\displaystyle\\left\\\{\\mathbf\{g\}^\{\(i\)\}\_\{1\},\\mathbf\{g\}^\{\(i\)\}\_\{2\},\\dots\\right\\\},\\qquad\\mathbf\{g\}^\{\(i\)\}\_\{\\ell\}=\\mathbf\{C\}\_\{i\}\\bigl\(\\mathbf\{A\}\_\{i\}\\bigr\)^\{\\ell\-1\}\\mathbf\{B\}\_\{i\}\.When an explicit delayτi\\tau\_\{i\}exists in the system, the Markov parameters exhibit a characteristic pattern: in the ideal noise\-free case, the firstτi\\tau\_\{i\}Markov parameters vanish, and nonzero responses begin at delayτi\+1\\tau\_\{i\}\+1\. This relationship between the Markov parameters and the time delay suggests that simply estimating the Markov parameters, without explicitly accounting for the delay, yields an estimate of the time\-delay system\. Furthermore, estimating the regimes represented by the time\-delay system can be reformulated as the problem of estimating mixtures of Markov parameters from the data\. Here, inLDSmeetsCTR, it was theoretically shown that mixtures of Markov parameters in observed data can be extracted by decomposing a properly generated third\-order tensor \(i\.e\., the system tensor\) as a superposition of rank\-1 tensors\. In this work, we construct an algorithm that applies this theoretical framework to regime estimation for time\-delay systems in real\-world streaming settings\. In the rest of this section, we focus on the proposed algorithm\.

![Refer to caption](https://arxiv.org/html/2605.26191v1/x1.png)Figure 1:Overview ofDelayMix\. Given the current data window𝐗c\\mathbf\{X\}^\{c\}and the current exogenous variables𝐔c\\mathbf\{U\}^\{c\}, the proposed method incrementally constructs a moment\-based system tensor, decomposes it into regime\-specific Markov parameters, realizes for each regime an equivalent delay\-free state\-space model, and uses these models to forecast future trajectories\.
### 2\.2Online Forecasting with Moment\-Based System Tensors

We now explain howDelayMixuses moment\-based system tensors for online modeling and the main idea is to keep an up\-to\-date empirical estimate of a system tensor at each update that summarizes the higher\-order moments of\(𝐱,𝐮\)\(\\mathbf\{x\},\\mathbf\{u\}\), and to periodically factorize this tensor to update a mixture of delay\-free state\-space models\. We then use these models to infer the current dynamical regime and to forecast future outputs based on the exogenous inputs\.

#### 2\.2\.1Data Windows and Exogenous Variables

We assume that the data stream𝐗\\mathbf\{X\}, denoted as\{𝐱​\(t\)\}t≥1\\\{\\mathbf\{x\}\(t\)\\\}\_\{t\\geq 1\}, arrives sequentially\. The algorithm runs at set update times, processing a sliding window of the most recent data\.

###### Definition 1\(Current data window:𝐗c\\mathbf\{X\}^\{c\}\)\.

Let𝐗c=𝐗\[t:tn\]\\mathbf\{X\}^\{c\}=\\mathbf\{X\}\[t:t\_\{n\}\]denote the subsequence of lengthlcl\_\{c\}extracted from the observed time series, starting at time indexttand ending attnt\_\{n\}\.

###### Definition 2\(Current exogenous variables:𝐔c\\mathbf\{U\}^\{c\}\)\.

Let𝐔c=𝐔\[t:tn\+ls\]\\mathbf\{U\}^\{c\}=\\mathbf\{U\}\[t:t\_\{n\}\+l\_\{s\}\]denote the subsequence of exogenous inputs available at the current update, wherelsl\_\{s\}is the forecasting horizon\. The segment𝐔\[t:tn\]\\mathbf\{U\}\[t:t\_\{n\}\]is used for updating the model, and𝐔\[tn\+1:tn\+ls\]\\mathbf\{U\}\[t\_\{n\}\+1:t\_\{n\}\+l\_\{s\}\]is used for forecasting\.

We focus on modeling streaming data as a mixture of time\-delay systems and on using dynamics for forecasting, which is one of the primary applications ofDelayMix\. At each update, the algorithm observes the current data window𝐗c\\mathbf\{X\}^\{c\}and the corresponding exogenous input sequence𝐔c\\mathbf\{U\}^\{c\}, along with the model estimated up to the previous update\. With this information, our goal is to refine the current system mixture and predict thelsl\_\{s\}\-step\-ahead outputs\.

We formalize the problem as follows\.

###### Problem 1\.

Giventhe current data window𝐗c\\mathbf\{X\}^\{c\}, the exogenous input sequence𝐔c\\mathbf\{U\}^\{c\}sampled at regular time intervals, the previously estimated model parametersθp\\theta\_\{p\}, and the previous system tensor𝒮p\\mathcal\{S\}\_\{p\},at each updatedo the following:

- •Update the system tensor𝒮\\mathcal\{S\}and the current system parameters θ=\{𝐀c,𝐁c,𝐂c\}\\displaystyle\\theta=\\bigl\\\{\\mathbf\{A\}\_\{c\},\\mathbf\{B\}\_\{c\},\\mathbf\{C\}\_\{c\}\\bigr\\\}to better fit the newly observed data;
- •Predict the future outputs𝐗f=𝐗\[tn\+1:tn\+ls\]\\mathbf\{X\}^\{f\}=\\mathbf\{X\}\[t\_\{n\}\+1:t\_\{n\}\+l\_\{s\}\]given the future inputs𝐔\[tn\+1:tn\+ls\]\\mathbf\{U\}\[t\_\{n\}\+1:t\_\{n\}\+l\_\{s\}\];
- •Maintain computational efficiency and bounded memory usage as new data arrives\.

To solve Problem[1](https://arxiv.org/html/2605.26191#Thmproblem1),DelayMixorganizes the online estimation and forecasting pipeline into three parts: DynamicMomentCollection \(system tensor collection\), ModelDatabaseAdaption \(online model adaptation via tensor factorization and system realization\), and FuturePrediction \(forecasting trajectories with the delay\-free models\)\. We describe each part in turn\. Figure[1](https://arxiv.org/html/2605.26191#S2.F1)shows an overview ofDelayMix, from the current data window and exogenous inputs to the updated models and forecasts\.

The goal of DynamicMomentCollection is to maintain an up\-to\-date empirical estimate of the global system tensor𝒮\\mathcal\{S\}using streaming input\-output data\. We build structured higher\-order moments to capture the temporal dependencies from the underlying switching dynamical systems\. Our approach uses a compact method inspired byLDSmeetsCTRthat does not require storing raw data\. DynamicMomentCollection keeps running estimates of the necessary moment blocks and updates𝒮\\mathcal\{S\}incrementally as new samples arrive\. This keeps memory usage low while preserving a tensor whose CP decomposition can recover the regime\-specific Markov parameters, up to scaling and permutation\. Concretely, for a sub\-window start timeτ\\tauand a triplet\(k1,k2,k3\)∈\{1,…,kmax\}3\(k\_\{1\},k\_\{2\},k\_\{3\}\)\\in\\\{1,\\dots,k\_\{\\max\}\\\}^\{3\}, define

t1=τ\+k1,t2=τ\+k1\+k2\+1,t3=τ\+k1\+k2\+k3\+2,t~1=τ,t~2=τ\+k1\+1,t~3=τ\+k1\+k2\+2\.t\_\{1\}=\\tau\+k\_\{1\},t\_\{2\}=\\tau\+k\_\{1\}\+k\_\{2\}\+1,t\_\{3\}=\\tau\+k\_\{1\}\+k\_\{2\}\+k\_\{3\}\+2,\\tilde\{t\}\_\{1\}=\\tau,\\tilde\{t\}\_\{2\}=\\tau\+k\_\{1\}\+1,\\tilde\{t\}\_\{3\}=\\tau\+k\_\{1\}\+k\_\{2\}\+2\.Using these indices, we form empirical sixth\-order moments over three output–input pairs and, for computational efficiency, we group each pair into a single mode usingvec​\(𝐱​\(tj\)​𝐮​\(tj\)⊤\)​\(j=1,2,3\)\\mathrm\{vec\}\(\\mathbf\{x\}\(t\_\{j\}\)\\mathbf\{u\}\(t\_\{j\}\)^\{\\top\}\)\(j=1,2,3\), yielding a third\-order moment blockℳ′​\(k1,k2,k3\)∈ℝp×p×p\\mathcal\{M\}^\{\\prime\}\(k\_\{1\},k\_\{2\},k\_\{3\}\)\\in\\mathbb\{R\}^\{p\\times p\\times p\}\. We combine these blocks across validτ\\tauand triplets to form the system tensor𝒮\\mathcal\{S\}using an incremental update\.

#### 2\.2\.2ModelDatabaseAdaption: Online Model Adaptation via Tensor Factorization

The goal of ModelDatabaseAdaption is to find and update the set of dynamical systemsθ\\thetafrom a moment\-based tensor representation\. While DynamicMomentCollection keeps updating𝒮\\mathcal\{S\}, extracting separate system modes requires a more involved tensor decomposition and system realization step\. This step is triggered only when a change in dynamics is detected or when enough new data have been collected\.

The adaptation procedure has the following steps:

1. 1\.Tensor decomposition:We apply the Alternating Least Squares \(ALS\) algorithm to the system tensor𝒮\\mathcal\{S\}to approximate it by a rank\-RRCP decomposition 𝒮≈∑i=1R𝐪i\(1\)⊗𝐪i\(2\)⊗𝐪i\(3\),\\displaystyle\\mathcal\{S\}\\approx\\sum\_\{i=1\}^\{R\}\\mathbf\{q\}^\{\(1\)\}\_\{i\}\\otimes\\mathbf\{q\}^\{\(2\)\}\_\{i\}\\otimes\\mathbf\{q\}^\{\(3\)\}\_\{i\},where\(𝐪i\(1\),𝐪i\(2\),𝐪i\(3\)\)\(\\mathbf\{q\}^\{\(1\)\}\_\{i\},\\mathbf\{q\}^\{\(2\)\}\_\{i\},\\mathbf\{q\}^\{\(3\)\}\_\{i\}\)are factor vectors associated with componentii\. Under the assumptions stated in Theorem[2](https://arxiv.org/html/2605.26191#Thmtheorem2), the rank\-one components can be associated with individual regimes and their stacked Markov parameters\.
2. 2\.Markov parameter reconstruction and delay\-free realization:For each componentii, we rearrange the factor vectors\(𝐪i\(1\),𝐪i\(2\),𝐪i\(3\)\)\(\\mathbf\{q\}^\{\(1\)\}\_\{i\},\\mathbf\{q\}^\{\(2\)\}\_\{i\},\\mathbf\{q\}^\{\(3\)\}\_\{i\}\)into a sequence of estimated Markov parameters\{𝐠^k\(i\)\}k=1K\\\{\\hat\{\\mathbf\{g\}\}^\{\(i\)\}\_\{k\}\\\}\_\{k=1\}^\{K\}by following the construction of the system tensor \(see also Section[3](https://arxiv.org/html/2605.26191#S3)\)\. The initial zero pattern in this sequence implicitly encodes the effective delay for regimeii, but we do not explicitly estimate the delay\. Instead, we treat\{𝐠^k\(i\)\}\\\{\\hat\{\\mathbf\{g\}\}^\{\(i\)\}\_\{k\}\\\}as the Markov parameters of an input–output equivalent system and apply a standard realization procedure, such as the Ho–Kalman algorithmHoKalmanto obtain a delay\-free state\-space model: \{𝐀i,𝐁i,𝐂i\}=Ho\-Kalman\(\{𝐠^k\(i\)\}k=1K\)\.\\displaystyle\\bigl\\\{\\mathbf\{A\}\_\{i\},\\mathbf\{B\}\_\{i\},\\mathbf\{C\}\_\{i\}\\bigl\\\}=\\text\{Ho\-Kalman\}\\bigl\(\\\{\\hat\{\\mathbf\{g\}\}^\{\(i\)\}\_\{k\}\\\}\_\{k=1\}^\{K\}\\bigr\)\.By Theorem[1](https://arxiv.org/html/2605.26191#Thmtheorem1)in Section[3](https://arxiv.org/html/2605.26191#S3), the resulting delay\-free model\{𝐀i,𝐁i,𝐂i\}\\bigl\\\{\\mathbf\{A\}\_\{i\},\\mathbf\{B\}\_\{i\},\\mathbf\{C\}\_\{i\}\\bigl\\\}is input–output equivalent to a time\-delay system and thus it faithfully captures the behavior of regimeii\. Finally, these parameters can be optionally fine\-tuned so as to better fit the current data𝐗c\\mathbf\{X\}^\{c\}, and the collection of all such systems constitutes the updated setθ\\theta\.

In summary, ModelDatabaseAdaption uses tensor factorization to separate regime\-specific Markov parameters from the system tensor and then realizes, for each regime, a delay\-free state\-space model that is suitable for inference and forecasting\.

#### 2\.2\.3FuturePrediction: Forecasting Trajectories with Delay\-Free Models

The goal of FuturePrediction is to infer a suitable initial state for the currently active system and use it, together with the learned models, to generate future values\. Given the updated set of delay\-free systemsθ\\thetaand the current data window\(𝐗c,𝐔c\)\(\\mathbf\{X\}^\{c\},\\mathbf\{U\}^\{c\}\), we first infer the active regime and its latent states, for example, by running a bank of Kalman filters and smoothers, one per regime\. We utilize the forward\-pass \(filter\) and backward\-pass \(smoother\) equations\.

The backward pass gives a smoothed estimate of the state trajectory, including an estimate of the current \(or initial\) state𝐳^​\(0\)\\hat\{\\mathbf\{z\}\}\(0\)for the active delay\-free system\. Using this state estimate and the identified system matrices\{𝐀i,𝐁i,𝐂i\}\\bigl\\\{\\mathbf\{A\}\_\{i\},\\mathbf\{B\}\_\{i\},\\mathbf\{C\}\_\{i\}\\bigl\\\}, we then simulate the delay\-free dynamics forward in time, using the future inputs𝐔\[tn\+1:tn\+ls\]\\mathbf\{U\}\[t\_\{n\}\+1:t\_\{n\}\+l\_\{s\}\], to generate thelsl\_\{s\}\-step\-ahead future values𝐗f\\mathbf\{X\}^\{f\}\.

By Theorem[1](https://arxiv.org/html/2605.26191#Thmtheorem1), each delay\-free model\{𝐀i,𝐁i,𝐂i\}\\bigl\\\{\\mathbf\{A\}\_\{i\},\\mathbf\{B\}\_\{i\},\\mathbf\{C\}\_\{i\}\\bigl\\\}obtained in this way is input–output equivalent to a time\-delay system representing regimeii\. This means forecasting with these delay\-free models is theoretically justified\.

## 3Theoretical Analysis

In this section, we explain the theoretical foundations behind the modeling approach\. We also discuss the assumptions and algorithmic design choices introduced earlier\. We focus on two main components that are central toDelayMix\. First, we represent each time\-delay MIMO system, which corresponds to one regime in the mixture model, using a delay\-free state\-space model built from its Markov parameters\. Second, we recover the regime\-specific Markov parameters from a moment\-based system tensor using CP decomposition\. We then analyze the computational complexity of the online algorithmDelayMix\.

### 3\.1Equivalence of a Time\-Delay System to a Standard MIMO Model

We begin by formalizing the relationship between a time\-delay MIMO system and an equivalent delay\-free state\-space model\. Recall from Section[2\.1](https://arxiv.org/html/2605.26191#S2.SS1)that each regimeiiin the mixture model is described by a time\-delay MIMO system\. For clarity, we first focus on a single system and drop the regime indexii\. The equivalence result in this subsection supports the realization step in ModelDatabaseAdaption and the forecasting step FuturePrediction inDelayMix\.

Consider a single segment of a time\-delay system described by

𝐳​\(t\+1\)\\displaystyle\\mathbf\{z\}\(t\+1\)=𝐀′​𝐳​\(t\)\+𝐁′​𝐮​\(t−τ\),\\displaystyle=\\mathbf\{A\}^\{\\prime\}\\mathbf\{z\}\(t\)\+\\mathbf\{B\}^\{\\prime\}\\mathbf\{u\}\(t\-\\tau\),\(3\)𝐱​\(t\)\\displaystyle\\mathbf\{x\}\(t\)=𝐂′​𝐳​\(t\),\\displaystyle=\\mathbf\{C\}^\{\\prime\}\\mathbf\{z\}\(t\),\(4\)whereτ≥0\\tau\\geq 0represents an explicit delay in the state update,𝐀′∈ℝk×k\\mathbf\{A\}^\{\\prime\}\\in\\mathbb\{R\}^\{k\\times k\}governs the latent dynamics,𝐁′∈ℝk×dc\\mathbf\{B\}^\{\\prime\}\\in\\mathbb\{R\}^\{k\\times d\_\{c\}\}captures the influence of the delayed input𝐮​\(t−τ\)\\mathbf\{u\}\(t\-\\tau\), and𝐂′∈ℝd×k\\mathbf\{C\}^\{\\prime\}\\in\\mathbb\{R\}^\{d\\times k\}maps the latent state to the observation\. For simplicity, we do not include direct\-feedthrough terms or additive noise\.

The Markov parameters of this delayed system are defined by its impulse response to an input that is zero everywhere except at a single time step\. They are given by

hj=\{0,for​1≤j≤τ,𝐂′​\(𝐀′\)j−τ−1​𝐁′,for​j\>τ\.\\displaystyle h\_\{j\}=\\begin\{cases\}0,&\\text\{for \}1\\leq j\\leq\\tau,\\\\ \\mathbf\{C\}^\{\\prime\}\(\\mathbf\{A\}^\{\\prime\}\)^\{j\-\\tau\-1\}\\mathbf\{B\}^\{\\prime\},&\\text\{for \}j\>\\tau\.\\end\{cases\}The firstτ\\tauzeros represent the delay, while the subsequent nonzero entries characterize the response after the delay\. It is well known that there exists a delay\-free state\-space realization whose Markov parameters match this sequence\. In other words, the explicit delay can be included in an augmented state\.

The next theorem shows that any minimal delay\-free realization with matching Markov parameters is input\-output equivalent to the original delayed system\.

###### Theorem 1\.

Suppose that the delayed system \([3](https://arxiv.org/html/2605.26191#S3.E3)\)–\([4](https://arxiv.org/html/2605.26191#S3.E4)\) is minimal \(reachable and observable\)\. Let\(𝐀,𝐁,𝐂\)\(\\mathbf\{A\},\\mathbf\{B\},\\mathbf\{C\}\)be any minimal delay\-free state\-space realization whose Markov parameters\{𝐂𝐀j−1​𝐁\}j≥1\\\{\\mathbf\{C\}\\mathbf\{A\}^\{j\-1\}\\mathbf\{B\}\\\}\_\{j\\geq 1\}coincide with those of the delayed system\. Then the delayed and delay\-free systems are input–output equivalent: for any input sequence and any initial state of the delayed system, there exists a corresponding initial state of the delay\-free realization such that the resulting output sequences coincide\.

Proof sketch of Theorem[1](https://arxiv.org/html/2605.26191#Thmtheorem1)\.The claim follows from the standard augmented\-state construction\. Define an augmented state that contains the original latent state and a buffer of the pastτ\\tauinputs\. The augmented transition matrix updates the original state using the oldest buffered input, shifts the buffer by one step, and inserts the current input into the buffer\. This gives a delay\-free state\-space realization\. By construction, its firstτ\\tauMarkov parameters are zero and the remaining Markov parameters coincide with those of the delayed system\. Therefore, for any input sequence, the two systems have the same input–output behavior under corresponding initial states\. ∎

Theorem[1](https://arxiv.org/html/2605.26191#Thmtheorem1)simplifies the problem of identifying a delayed system to estimating its possibly sparse sequence of Markov parameters from data\. This supports working in the standard state\-space framework without explicitly modeling delays, as long as the Markov parameters are correctly captured\. In particular, this supports the modeling choice inDelayMixto perform inference forecasting entirely in terms of delay\-free state\-space models realized from estimated Markov parameters\.

### 3\.2Moment Tensor Representation and Tensor Factorization

Next, we summarize the moment\-based tensor representation used in the online estimation algorithm ofDelayMix\. In the previous section, we introduced the system tensor𝒮\\mathcal\{S\}as a compact way to encode higher\-order moments of the input\-output data and as the object for CP decomposition in DynamicMomentCollection and ModelDatabaseAdaption\. Here we state a concise result fromLDSmeetsCTR, which shows that factorizing a properly constructed moment tensor is enough to recover the Markov parameters of the underlying systems\.

###### Theorem 2\(Recovery of Markov parameters via moment tensor factorizationLDSmeetsCTR\)\.

Let\(𝐮​\(t\),𝐱​\(t\)\)t≥1\(\\mathbf\{u\}\(t\),\\mathbf\{x\}\(t\)\)\_\{t\\geq 1\}be generated by a \(mixture of\) linear time\-delay MIMO systems with regime\-specific Markov parameters\{hj\(i\)\}j≥1\\\{h^\{\(i\)\}\_\{j\}\\\}\_\{j\\geq 1\}\. Assume that the exogenous inputs are independent and identically distributed, with a distribution that satisfies suitable moment conditions, and that the resulting moment tensor is well\-defined\. Then there exists a polynomial\-time procedure that constructs a tensor𝒮\\mathcal\{S\}from empirical moments of\(𝐮​\(t\),𝐱​\(t\)\)\(\\mathbf\{u\}\(t\),\\mathbf\{x\}\(t\)\)such that, under mild identifiability conditions, any sufficiently accurate CP decomposition of𝒮\\mathcal\{S\}recovers the Markov parameters\{hj\(i\)\}j≥1\\\{h^\{\(i\)\}\_\{j\}\\\}\_\{j\\geq 1\}of the constituent systems, up to permutation and scaling of the components\.

###### Proof\.

SeeLDSmeetsCTRfor the complete proof\. ∎

This theorem gives the theoretical justification for the tensor factorization step inDelayMix\. The CP factors obtained by applying ALS to the empirical system tensor in ModelDatabaseAdaption can be seen as estimates of the regime\-specific Markov parameters\. Combined with Theorem[1](https://arxiv.org/html/2605.26191#Thmtheorem1), this ensures that the subsequent realization step produces delay\-free models whose input\-output behavior matches that of the underlying time\-delay systems\.

![Refer to caption](https://arxiv.org/html/2605.26191v1/x2.png)Figure 2:Forecasting performance comparison\. On the delayed mixture benchmark,DelayMixoutperforms competing methods in terms of mean squared error \(MSE\) and mean absolute error \(MAE\)\.
### 3\.3Computational Complexity ofDelayMix

Finally, we analyze the computational complexity of the proposed online algorithmDelayMix\. The dominant computational cost arises from the ALS\-based tensor decomposition used to update the system tensor representation, along with the following realization of state\-space models\.

###### Lemma1\(Time complexity ofDelayMix\)\.

Given a new incoming tensor \(i\.e\., an updated system tensor𝒮\\mathcal\{S\}\), the time complexity of one update ofDelayMixis

𝒪​\(i​\(8​s3​d3​dc3​R\+6​s​d​dc​R2\)\+k3​lc\+k2​dc​ls\+k2​ls\+k​ls​d\),\\displaystyle\\mathcal\{O\}\\bigl\(i\\,\(8s^\{3\}d^\{3\}d\_\{c\}^\{3\}R\+6sdd\_\{c\}\\,R^\{2\}\)\+k^\{3\}l\_\{c\}\+k^\{2\}d\_\{c\}l\_\{s\}\+k^\{2\}l\_\{s\}\+kl\_\{s\}d\\bigr\),whereiiis the number of ALS iterations for tensor decomposition,ssis the maximum lag used to construct𝒮\\mathcal\{S\},RRis the CP rank,lcl\_\{c\}is the length of the current data window, andlsl\_\{s\}is the forecasting horizon\.

Proof sketch of Lemma[1](https://arxiv.org/html/2605.26191#Thmlemma1)\.Letp=2​s​d​dcp=2sdd\_\{c\}be the mode size of the constructed third\-order system tensor\. One ALS iteration for a rank\-RRCP decomposition costsO​\(p3​R\+p​R2\)O\(p^\{3\}R\+pR^\{2\}\), and repeating it foriiiterations gives the first term\. The remaining terms are the costs of state inference on the current window and simulation over the forecasting horizon\. ∎

The computational cost ofDelayMixmainly depends on the ALS\-based tensor decomposition, especially on the number of iterationsiineeded for convergence\. To reduce this cost in practice, the full decomposition step is only triggered when there is a significant mismatch between the current model and the newly updated system tensor is detected\. This helps avoid unnecessary updates\. Also, because the data changes gradually over time, consecutive system tensors are similar\. As a result, initializing ALS with the previous factor estimates usually leads to fast convergence \(smallii\)\. These properties allowDelayMixto work efficiently in streaming settings while maintaining accurate modeling and forecasting performance\.

## 4Experimental Results

In this section, we describe the performance ofDelayMixusing real datasets\. The experiments were designed to answer the following questions aboutDelayMix:

Q1\.Accuracy:How accurately does it predict future events?

Q2\.Scalability:How does it scale in terms of computational time and memory?

### 4\.1Experimental Setup

![Refer to caption](https://arxiv.org/html/2605.26191v1/x3.png)Figure 3:Critical difference diagram of datasets in terms of MSE and MAE\.Datasets\.To evaluateDelayMix, we use three real\-world datasets\.

For the evaluation of predictive accuracy \(Q1\) and scalability \(Q2\), we used the following three publicly available real\-world datasets:

- •Ship\-OOD/Ship\-INDShipdatasetis a 4\-DOF ship maneuvering dataset under environmental disturbances\. The 4 exogenous variables include propeller speed, rudder angles, and wind; outputs are also 4 \(linear and angular velocities\)\. We use both in\-distribution \(Ship\-IND\) and out\-of\-distribution \(Ship\-OOD\) test sets\.
- •RobotRobotdatasetis a dataset of a 6\-DOF industrial robot arm, with 6 exogenous inputs \(motor torques\) and 6 outputs \(joint angles\)\.

Baseline Methods\.We compareDelayMixwith state\-of\-the\-art models that explicitly incorporate exogenous variables, including StableReLiNet, ReLiNetReLiNet, and TimeXerTimeXer\. Many existing forecasting models either omit exogenous inputs or treat them as endogenous, limiting their ability to model structured input\-output relationships\. Accordingly, we exclude such models from our main comparison\.

Hyperparameter Setting ofDelayMix\.We selectedρ∈\{0\.5,0\.6,0\.7,0\.8,0\.9,1\.0\}\\rho\\in\\\{0\.5,0\.6,0\.7,0\.8,0\.9,1\.0\\\}andR∈\{2,4,8,10\}R\\in\\\{2,4,8,10\\\}using the validation set, and fixeds=3s=3\.

### 4\.2Q1: Accuracy

We evaluatedDelayMixusing mean squared error \(MSE\) and mean absolute error \(MAE\), where lower values indicate better accuracy\. Results are averaged over five trials for robustness\. As shown in Fig\.[2](https://arxiv.org/html/2605.26191#S3.F2),DelayMixconsistently yields lower errors than baselines\. The critical difference diagram in Fig\.[3](https://arxiv.org/html/2605.26191#S4.F3), based on the Wilcoxon\-Holm testWilcoxonHolmmethod, confirms these improvements are statistically significant\. Table[1](https://arxiv.org/html/2605.26191#S4.T1)reports how often each method ranks best, further highlighting the advantage ofDelayMix\.

A key strength ofDelayMixis its ability to adaptively estimate a mixture of time\-delay systems in streaming settings, enabling it to capture dynamic and heterogeneous temporal patterns without prior knowledge of delay structures or access to past data\. While TimeXer performs well on the stationary Ship\-IND dataset, it falls short on the non\-stationary Ship\-OOD data, where adaptive modeling is crucial\. In contrast, models like StableReLiNet and ReLiNet tend to degrade under streaming updates, likely due to unstable parameter adaptation\. For example, ReLiNet shows noticeably lower accuracy in several scenarios \(see Fig\.[2](https://arxiv.org/html/2605.26191#S3.F2)\)\. These results show thatDelayMixenables accurate real\-time prediction in dynamic streaming settings\.

Table 1:1st1^\{\\text\{st\}\}Count in dataset \(higher is better\)\.
### 4\.3Q2: Scalability

Next, we evaluated the performance ofDelayMixin terms of computation time\. The left column of Fig\.[4](https://arxiv.org/html/2605.26191#S4.F4)shows the wall clock time of the experiments performed on each dataset\. Thanks to an effective parameter estimation algorithm, the computation time ofDelayMixis independent of the length of the data stream\. Several small spikes are due to the execution of ModelDatabaseAdaption\. CP decomposition for updating the regime database increases the computation time for that time step\. The right column of Fig\.[4](https://arxiv.org/html/2605.26191#S4.F4)shows the average computation time for the entire tensor stream\. Notably, TimeXer and ReLiNet were run in a relatively fast environment, as all computations, including forecasting future values and parameter updating using lightweight single\-sample gradients, are performed on the GPU\. However,DelayMixis still overwhelmingly faster than them\. This demonstrates thatDelayMixachieves a high trade\-off between accuracy and efficiency\.

![Refer to caption](https://arxiv.org/html/2605.26191v1/x4.png)Figure 4:Efficiency of streamingDelayMix: Our method is consistently faster than other baseline methods\.

## 5Related Work

This section reviews related work on mixtures of dynamical systems and time\-series forecasting\.

Learning Mixtures of Dynamics\.Traditionally, dynamical system estimation has centered on a single systemSINDy;ode1;si\_icml\_18;KDD2023\_PDE;mio\-sindy;LaNoLem;ICLR2025\_si;ICMR2025\_si, but recent advances enable learning mixed systems—where multiple dynamical regimes coexist—from unlabeled trajectory data\. For example,MixtureofLDSproposes a spectral algorithm that learns a mixture from short unlabeled trajectories with end\-to\-end guarantees, andLDSmeetsCTRintroduces a tensor\-based framework that recovers LDS mixtures via Markov parameters estimated from high\-order moments\. These approaches are theoretically appealing but inherently offline: constructing the required moment tensors assumes access to the full dataset, which is impractical in streaming scenarios such as industrial process control\. Moreover, they handle delays only implicitly and do not explicitly model unknown or time\-varying delays, making it difficult to adapt to abrupt changes in system dynamics\.DelayMixbuilds on the framework ofLDSmeetsCTRand extends it to streaming data by incrementally updating a system tensor and adaptively estimating multiple time\-delay systems from incoming observations\.

Time\-Series Forecasting Methods\.Classical forecasting methods such as ARIMA and Kalman filtersARIMA;LDS;Dynammoremain standard tools, and recent deep learning modelsInformer;DeepAR;DLinear;PatchTST;Fredformerhave achieved strong empirical performance\. However, many of these methods either ignore exogenous variables or treat them in the same way as endogenous variables, although exogenous inputs often play a distinct and crucial role in prediction\. To better exploit exogenous information, ARIMAXARIMAXextends ARIMA with exogenous regressors, while recent neural approaches such as ReLiNetReLiNetand TimeXerTimeXerincorporate exogenous signals through switched linear systems and Transformer\-based attention, respectively, achieving competitive results\. However, these models generally do not model time delays explicitly and tend to struggle in streaming settings due to high computational cost and frequent retraining, whereasDelayMixis designed for online operation: it maintains a compact summary of input–output relationships, explicitly captures delay dynamics, and updates models efficiently in real time\.

## 6Conclusion

We introduceDelayMix, an online framework for modeling dynamic mixtures of time\-delay systems from streaming time series\. Our approach addresses key challenges in streaming scenarios, including nonstationarity, unknown delays, and limited resources, while maintaining both accuracy and scalability over long horizons\.DelayMixincrementally builds and decomposes high\-order moment tensors to extract regime\-specific Markov parameter sequences, which efficiently capture system dynamics and input\-output delays in a compact, delay\-sensitive representation\. This allows for delay\-aware identification even when regimes switch, without assuming known change points or fixed delay structures\. The framework does not require direct access to the system or prior knowledge of delay lengths\. Instead, it uses moment\-based summaries that are updated online, so it does not need to store the entire history, keeping memory use low and computation efficient\.DelayMixcan also update its set of candidate systems when the current model no longer fits new observations, helping it adapt to changing dynamics while avoiding unbounded memory growth or switching costs\. Our experiments with real\-world datasets show thatDelayMixconsistently achieves better prediction accuracy and computational efficiency than other methods by accurately tracking changing delays and adapting to regime shifts\.

## Acknowledgements

This work was supported by JSPS KAKENHI Grant\-in\-Aid for Scientific Research Number JP26H02499, JST CREST JPMJCR23M3, JST K Program JPMJKP25Y6, JST COI\-NEXT JPMJPF2009, JST COI\-NEXT JPMJPF2115, Future Social Value Co\-Creation Project \- Osaka University \.

## References

## Appendix

## Appendix ASymbols and Definitions

Table 2:Symbols and definitions\.Table[2](https://arxiv.org/html/2605.26191#A1.T2)lists the symbols and definitions used in this paper\.

## Appendix BDetails of algorithms

### B\.1Ho\-Kalman algorithm

Here, we show the overrview of Ho\-Kalman algorithmHoKalman\.

Algorithm 1Ho\-kalman \(𝐆\\mathbf\{G\}\)1:Input:\(a\) A markov parameter

𝐆=\{𝐠1,…,𝐠2​s\}\\mathbf\{G\}=\\\{\\mathbf\{g\}\_\{1\},\.\.\.,\\mathbf\{g\}\_\{2s\}\\\}
2:Output:System parameters

\{𝐀,𝐁,𝐂\}\\\{\\mathbf\{A\},\\mathbf\{B\},\\mathbf\{C\}\\\}
3:Form the Hankel matrix

H^∈ℝm​s×p​\(s\+1\)\\hat\{H\}\\in\\mathbb\{R\}^\{ms\\times p\(s\+1\)\}from

𝐆\\mathbf\{G\}as

H^=\[𝐠1𝐠2⋯𝐠s\+1𝐠2𝐠3⋯𝐠s\+2⋮⋮⋱⋮𝐠s𝐠s\+1⋯𝐠2​s\]\\hat\{H\}=\\begin\{bmatrix\}\\mathbf\{g\}\_\{1\}&\\mathbf\{g\}\_\{2\}&\\cdots&\\mathbf\{g\}\_\{s\+1\}\\\\ \\mathbf\{g\}\_\{2\}&\\mathbf\{g\}\_\{3\}&\\cdots&\\mathbf\{g\}\_\{s\+2\}\\\\ \\vdots&\\vdots&\\ddots&\\vdots\\\\ \\mathbf\{g\}\_\{s\}&\\mathbf\{g\}\_\{s\+1\}&\\cdots&\\mathbf\{g\}\_\{2s\}\\end\{bmatrix\}
4:

H^−∈ℝm​s×p​s←\\hat\{H\}^\{\-\}\\in\\mathbb\{R\}^\{ms\\times ps\}\\leftarrowfirst

p​spscolumns of

H^\\hat\{H\}
5:

L^∈ℝm​s×p​s←\\hat\{L\}\\in\\mathbb\{R\}^\{ms\\times ps\}\\leftarrowrank

nnapproximation of

H^−\\hat\{H\}^\{\-\}via SVD

6:

U,Σ,V=SVD​\(L^\)U,\\Sigma,V=\\mathrm\{SVD\}\(\\hat\{L\}\)
7:

O^∈ℝm​s×n←U​Σ1/2\\hat\{O\}\\in\\mathbb\{R\}^\{ms\\times n\}\\leftarrow U\\Sigma^\{1/2\}
8:

Q^∈ℝn×p​s←Σ1/2​V⊤\\hat\{Q\}\\in\\mathbb\{R\}^\{n\\times ps\}\\leftarrow\\Sigma^\{1/2\}V^\{\\top\}
9:

C^←\\hat\{C\}\\leftarrowfirst

mmrows of

O^\\hat\{O\}
10:

B^←\\hat\{B\}\\leftarrowfirst

ppcolumns of

Q^\\hat\{Q\}
11:

H^\+∈ℝm​s×p​s←\\hat\{H\}^\{\+\}\\in\\mathbb\{R\}^\{ms\\times ps\}\\leftarrowlast

p​spscolumns of

H^\\hat\{H\}
12:

A^←O^†​H^\+​Q^†\\hat\{A\}\\leftarrow\\hat\{O\}^\{\\dagger\}\\hat\{H\}^\{\+\}\\hat\{Q\}^\{\\dagger\}
13:return

A^,B^,C^\\hat\{A\},\\hat\{B\},\\hat\{C\}

### B\.2Algorithm Overview

Here, we show the overrview ofDelayMix\.

Algorithm 2DelayMix\(𝐗c,𝐔c,𝒮p,θp\\mathbf\{X\}^\{c\},\\mathbf\{U\}^\{c\},\\mathcal\{S\}\_\{p\},\\theta\_\{p\}\)1:Input:\(a\) Current data window

𝐗c=\{𝐱​\(t\)\}\\mathbf\{X\}^\{c\}=\\\{\\mathbf\{x\}\(t\)\\\}\(b\) Current exogenous window

𝐔c=\{𝐮​\(t\)\}\\mathbf\{U\}^\{c\}=\\\{\\mathbf\{u\}\(t\)\\\}\(c\) Previous system tensor

𝒮p\\mathcal\{S\}\_\{p\}\(d\) Previous system parameters

θp\\theta\_\{p\}
2:Output:\(a\) Current system parameters

θ\\theta\(b\) Updated system tensor

𝒮\\mathcal\{S\}\(c\) Future values

𝐗f\\mathbf\{X\}^\{f\}
3:/\* DynamicMomentCollection: update global system tensor \*/

4:

𝒮←𝒮p\\mathcal\{S\}\\leftarrow\\mathcal\{S\}\_\{p\}
5:Set

kmaxk\_\{\\max\}, and let

p=d​dcp=dd\_\{c\},

D=kmax​pD=k\_\{\\max\}p
6:Let

𝒥​\(k\)\\mathcal\{J\}\(k\)denote the index range

\{\(k−1\)​p\+1,…,k​p\}\\\{\(k\-1\)p\+1,\\dots,kp\\\}
7:for all

\(k1,k2,k3\)∈\{1,…,kmax\}3\(k\_\{1\},k\_\{2\},k\_\{3\}\)\\in\\\{1,\\dots,k\_\{\\max\}\\\}^\{3\}do

8:for alladmissible sub\-window starts

τ\\tauin the current windowdo

9:

t1←τ\+k1t\_\{1\}\\leftarrow\\tau\+k\_\{1\}
10:

t2←τ\+k1\+k2\+1t\_\{2\}\\leftarrow\\tau\+k\_\{1\}\+k\_\{2\}\+1
11:

t3←τ\+k1\+k2\+k3\+2t\_\{3\}\\leftarrow\\tau\+k\_\{1\}\+k\_\{2\}\+k\_\{3\}\+2
12:

t~1←τ\\tilde\{t\}\_\{1\}\\leftarrow\\tau,

t~2←τ\+k1\+1\\tilde\{t\}\_\{2\}\\leftarrow\\tau\+k\_\{1\}\+1,

t~3←τ\+k1\+k2\+2\\tilde\{t\}\_\{3\}\\leftarrow\\tau\+k\_\{1\}\+k\_\{2\}\+2
13:

𝐦\(1\)←vec​\(𝐱​\(t1\)​𝐮​\(t~1\)⊤\)\\mathbf\{m\}^\{\(1\)\}\\leftarrow\\mathrm\{vec\}\\\!\\left\(\\mathbf\{x\}\(t\_\{1\}\)\\mathbf\{u\}\(\\tilde\{t\}\_\{1\}\)^\{\\top\}\\right\)
14:

𝐦\(2\)←vec​\(𝐱​\(t2\)​𝐮​\(t~2\)⊤\)\\mathbf\{m\}^\{\(2\)\}\\leftarrow\\mathrm\{vec\}\\\!\\left\(\\mathbf\{x\}\(t\_\{2\}\)\\mathbf\{u\}\(\\tilde\{t\}\_\{2\}\)^\{\\top\}\\right\)
15:

𝐦\(3\)←vec​\(𝐱​\(t3\)​𝐮​\(t~3\)⊤\)\\mathbf\{m\}^\{\(3\)\}\\leftarrow\\mathrm\{vec\}\\\!\\left\(\\mathbf\{x\}\(t\_\{3\}\)\\mathbf\{u\}\(\\tilde\{t\}\_\{3\}\)^\{\\top\}\\right\)
16:

𝒮​\[𝒥​\(k1\),𝒥​\(k2\),𝒥​\(k3\)\]←𝒮​\[𝒥​\(k1\),𝒥​\(k2\),𝒥​\(k3\)\]\+𝐦\(3\)⊗𝐦\(2\)⊗𝐦\(1\)\\mathcal\{S\}\[\\mathcal\{J\}\(k\_\{1\}\),\\mathcal\{J\}\(k\_\{2\}\),\\mathcal\{J\}\(k\_\{3\}\)\]\\leftarrow\\mathcal\{S\}\[\\mathcal\{J\}\(k\_\{1\}\),\\mathcal\{J\}\(k\_\{2\}\),\\mathcal\{J\}\(k\_\{3\}\)\]\+\\mathbf\{m\}^\{\(3\)\}\\otimes\\mathbf\{m\}^\{\(2\)\}\\otimes\\mathbf\{m\}^\{\(1\)\}
17:endfor

18:endfor

19:/\* Compute the error between the current system and the current data \*/

20:Compute predictions

\{𝐱^​\(j\)\}\\\{\\mathbf\{\\hat\{x\}\}\(j\)\\\}from

θp\\theta\_\{p\}on the current window \(e\.g\., via Kalman filtering; see Appendix[B\.3](https://arxiv.org/html/2605.26191#A2.SS3)\)\.

21:

Lc←∑j=tt\+tn‖𝐱​\(j\)−𝐱^​\(j\)‖2L\_\{c\}\\leftarrow\\sum^\{t\+t\_\{n\}\}\_\{j=t\}\\\|\\mathbf\{x\}\(j\)\-\\mathbf\{\\hat\{x\}\}\(j\)\\\|\_\{2\}
22:/\* The system

θp\\theta\_\{p\}is not updated if it sufficiently fits the current data\. \*/

23:if

Lc<ρL\_\{c\}<\\rhothen

24:

θ←θp\\theta\\leftarrow\\theta\_\{p\}
25:else

26:/\* ModelDatabaseAdaption \*/

27:

\{𝐪i\(1\),𝐪i\(2\),𝐪i\(3\)\}i=1R←\\\{\\mathbf\{q\}^\{\(1\)\}\_\{i\},\\mathbf\{q\}^\{\(2\)\}\_\{i\},\\mathbf\{q\}^\{\(3\)\}\_\{i\}\\\}\_\{i=1\}^\{R\}\\leftarrowCP\-Decomposition\(

𝒮\\mathcal\{S\}\)

28:for

i=1i=1to

RRdo

29:/\* Construct Markov parameters

𝐆i\\mathbf\{G\}\_\{i\}LDSmeetsCTR\*/

30:

𝒯←𝐪i\(1\)⊗𝐪i\(2\)⊗𝐪i\(3\)\\mathcal\{T\}\\leftarrow\\mathbf\{q\}^\{\(1\)\}\_\{i\}\\otimes\\mathbf\{q\}^\{\(2\)\}\_\{i\}\\otimes\\mathbf\{q\}^\{\(3\)\}\_\{i\}
31:Compute the Frobenius norm of

𝒯\\mathcal\{T\}along modes 2–3 and obtain

𝐯i\\mathbf\{v\}\_\{i\}
32:

𝐆i←\\mathbf\{G\}\_\{i\}\\leftarrowReshape

\(𝐯i/‖𝐯i‖2/3\)\\\!\\left\(\\mathbf\{v\}\_\{i\}/\\\|\\mathbf\{v\}\_\{i\}\\\|^\{2/3\}\\right\)
33:

\{𝐀i,𝐁i,𝐂i\}←\\\{\\mathbf\{A\}\_\{i\},\\mathbf\{B\}\_\{i\},\\mathbf\{C\}\_\{i\}\\\}\\leftarrowHo\-Kalman\(

𝐆i\\mathbf\{G\}\_\{i\}\)

34:endfor

35:/\* Identify the system that best explains

𝐗c\\mathbf\{X\}^\{c\}\*/

36:

θ←arg​min\{𝐀i,𝐁i,𝐂i\}​∑j=tt\+tn‖𝐱​\(j\)−𝐱^​\(j\)‖2\\theta\\leftarrow\\underset\{\\\{\\mathbf\{A\}\_\{i\},\\mathbf\{B\}\_\{i\},\\mathbf\{C\}\_\{i\}\\\}\}\{\\operatorname\*\{\\text\{arg\}\\,\\text\{min\}\}\}\\ \\sum^\{t\+t\_\{n\}\}\_\{j=t\}\\\|\\mathbf\{x\}\(j\)\-\\mathbf\{\\hat\{x\}\}\(j\)\\\|\_\{2\}
37:endif

38:/\* FuturePrediction \*/

39:Compute

𝐳^​\(0\)\\hat\{\\mathbf\{z\}\}\(0\)using

θ\\theta\(see Appendix[B\.3](https://arxiv.org/html/2605.26191#A2.SS3)\)\.

40:Generate future values

𝐗f\\mathbf\{X\}^\{f\}from

𝐳^​\(0\)\\hat\{\\mathbf\{z\}\}\(0\)and future inputs\.

41:return

θ,𝒮,𝐗f\\theta,\\mathcal\{S\},\\mathbf\{X\}^\{f\}

### B\.3Generate: Forecasting

Here, we describe the following the forward\-pass \(filter\) and backward\-pass \(smoother\) equations\.

##### Forward Pass \(Kalman Filter\):

𝝁^​\(t\)\\displaystyle\\hat\{\\boldsymbol\{\\mu\}\}\(t\)=𝐀​𝝁​\(t−1\)\+𝐁𝐮​\(t−1\),\\displaystyle=\\mathbf\{A\}\\boldsymbol\{\\mu\}\(t\-1\)\+\\mathbf\{B\}\\mathbf\{u\}\(t\-1\),\(5\)𝐏^​\(t\)\\displaystyle\\hat\{\\mathbf\{P\}\}\(t\)=𝐀𝐏​\(t−1\)​𝐀⊤\+𝚪,\\displaystyle=\\mathbf\{A\}\\mathbf\{P\}\(t\-1\)\\mathbf\{A\}^\{\\top\}\+\\mathbf\{\\Gamma\},\(6\)𝐊​\(t\)\\displaystyle\\mathbf\{K\}\(t\)=𝐏^​\(t\)​𝐂⊤​\(𝐂​𝐏^​\(t\)​𝐂⊤\+𝐑\)−1,\\displaystyle=\\hat\{\\mathbf\{P\}\}\(t\)\\mathbf\{C\}^\{\\top\}\(\\mathbf\{C\}\\hat\{\\mathbf\{P\}\}\(t\)\\mathbf\{C\}^\{\\top\}\+\\mathbf\{R\}\)^\{\-1\},\(7\)𝝁​\(t\)\\displaystyle\\boldsymbol\{\\mu\}\(t\)=𝝁^​\(t\)\+𝐊​\(t\)​\(𝐱​\(t\)−𝐂​𝝁^​\(t\)\),\\displaystyle=\\hat\{\\boldsymbol\{\\mu\}\}\(t\)\+\\mathbf\{K\}\(t\)\(\\mathbf\{x\}\(t\)\-\\mathbf\{C\}\\hat\{\\boldsymbol\{\\mu\}\}\(t\)\),\(8\)𝐏​\(t\)\\displaystyle\\mathbf\{P\}\(t\)=\(𝐈−𝐊​\(t\)​𝐂\)​𝐏^​\(t\)\.\\displaystyle=\(\\mathbf\{I\}\-\\mathbf\{K\}\(t\)\\mathbf\{C\}\)\\hat\{\\mathbf\{P\}\}\(t\)\.\(9\)

##### Backward Pass \(RTS Smoother\):

𝐕​\(t\)\\displaystyle\\mathbf\{V\}\(t\)=𝐏​\(t\)​𝐀⊤​𝐏^​\(t\+1\)−1,\\displaystyle=\\mathbf\{P\}\(t\)\\mathbf\{A\}^\{\\top\}\\hat\{\\mathbf\{P\}\}\(t\+1\)^\{\-1\},\(10\)𝐳^​\(t\)\\displaystyle\\hat\{\\mathbf\{z\}\}\(t\)=𝝁​\(t\)\+𝐕​\(t\)​\(𝐳^​\(t\+1\)−𝝁^​\(t\+1\)\)\.\\displaystyle=\\boldsymbol\{\\mu\}\(t\)\+\\mathbf\{V\}\(t\)\(\\hat\{\\mathbf\{z\}\}\(t\+1\)\-\\hat\{\\boldsymbol\{\\mu\}\}\(t\+1\)\)\.\(11\)

## Appendix CProofs

### C\.1Proof of Theorem[1](https://arxiv.org/html/2605.26191#Thmtheorem1)

###### Lemma2\(Delay\-free realization from Markov parameters\)\.

Let the Markov parameters of the delayed system \([3](https://arxiv.org/html/2605.26191#S3.E3)\)–\([4](https://arxiv.org/html/2605.26191#S3.E4)\) be given by\{gj\}j≥1\\\{g\_\{j\}\\\}\_\{j\\geq 1\}\. Then there exists a delay\-free state\-space model

𝐳~​\(t\+1\)\\displaystyle\\tilde\{\\mathbf\{z\}\}\(t\+1\)=𝐀​𝐳~​\(t\)\+𝐁𝐮​\(t\),\\displaystyle=\\mathbf\{A\}\\tilde\{\\mathbf\{z\}\}\(t\)\+\\mathbf\{B\}\\mathbf\{u\}\(t\),\(12\)𝐱^​\(t\)\\displaystyle\\mathbf\{\\hat\{x\}\}\(t\)=𝐂​𝐳~​\(t\),\\displaystyle=\\mathbf\{C\}\\tilde\{\\mathbf\{z\}\}\(t\),\(13\)with possibly higher\-dimensional latent state𝐳~​\(t\)\\tilde\{\\mathbf\{z\}\}\(t\)and suitable matrices𝐀,𝐁,𝐂\\mathbf\{A\},\\mathbf\{B\},\\mathbf\{C\}such that

gj=𝐂𝐀j−1​𝐁for all​j≥1\.\\displaystyle g\_\{j\}=\\mathbf\{C\}\\mathbf\{A\}^\{j\-1\}\\mathbf\{B\}\\quad\\text\{for all \}j\\geq 1\.In other words, the delay\-free model reproduces exactly the same Markov parameters as the original delayed system\.

###### Proof\.

Because the original system has an explicit input delay ofτ\\tausteps, the input𝐮​\(t\)\\mathbf\{u\}\(t\)only affects the output from timet\+τt\+\\tauonward\. Its input–output behavior is therefore fully characterized by the Markov parameter sequence\{gj\}j≥1\\\{g\_\{j\}\\\}\_\{j\\geq 1\}, which defines the impulse response

𝐱^​\(t\)=∑j=τ∞gj​𝐮​\(t−j\)\.\\displaystyle\\mathbf\{\\hat\{x\}\}\(t\)=\\sum\_\{j=\\tau\}^\{\\infty\}g\_\{j\}\\,\\mathbf\{u\}\(t\-j\)\.
We construct an equivalent delay\-free state\-space model by embedding the delayed inputs into an extended state vector\. Specifically, define the extended state to stack the lastτ\\tauinputs:

𝐳~​\(t\):=\[𝐮​\(t−τ\+1\)⊤𝐮​\(t−τ\+2\)⊤⋯𝐮​\(t\)⊤\]⊤\.\\displaystyle\\tilde\{\\mathbf\{z\}\}\(t\):=\\begin\{bmatrix\}\\mathbf\{u\}\(t\-\\tau\+1\)^\{\\top\}&\\mathbf\{u\}\(t\-\\tau\+2\)^\{\\top\}&\\cdots&\\mathbf\{u\}\(t\)^\{\\top\}\\end\{bmatrix\}^\{\\top\}\.With this construction, the delay is absorbed into the state update, and the output𝐱^​\(t\)\\mathbf\{\\hat\{x\}\}\(t\)can be written as a linear function of𝐳~​\(t\)\\tilde\{\\mathbf\{z\}\}\(t\)without explicit delay\.

As a result, the convolution representation induced by\{gj\}\\\{g\_\{j\}\\\}can be realized by a standard delay\-free state\-space model\. Classical realization techniques such as the Ho–Kalman algorithm or subspace methods applied to the Markov parameter sequence\{gj\}\\\{g\_\{j\}\\\}yield matrices\(𝐀,𝐁,𝐂\)\(\\mathbf\{A\},\\mathbf\{B\},\\mathbf\{C\}\)that satisfy

gj=𝐂𝐀j−1​𝐁for all​j≥1\.\\displaystyle g\_\{j\}=\\mathbf\{C\}\\mathbf\{A\}^\{j\-1\}\\mathbf\{B\}\\quad\\text\{for all \}j\\geq 1\.Hence the delayed system and the constructed delay\-free system share identical Markov parameters, which proves the claim\. ∎

###### Proof of Theorem[1](https://arxiv.org/html/2605.26191#Thmtheorem1)\.

By Lemma[2](https://arxiv.org/html/2605.26191#Thmlemma2), the delayed system \([3](https://arxiv.org/html/2605.26191#S3.E3)\)–\([4](https://arxiv.org/html/2605.26191#S3.E4)\) admits a delay\-free state\-space realization whose Markov parameters\{gj\}j≥1\\\{g\_\{j\}\\\}\_\{j\\geq 1\}coincide with those of the delayed system\. The input–output map of a linear time\-invariant system is completely determined by its Markov parameters \(equivalently, by its transfer function\)\. Therefore any two minimal state\-space realizations with the same Markov parameters are related by a similarity transformation and produce identical outputs for any input sequence\. Since both the delayed system and\(𝐀,𝐁,𝐂\)\(\\mathbf\{A\},\\mathbf\{B\},\\mathbf\{C\}\)are minimal and share the same Markov parameters, they are input–output equivalent\. ∎

### C\.2Proof of Lemma[1](https://arxiv.org/html/2605.26191#Thmlemma1)

The most time\-consuming part in DynamicMomentCollection and ModelDatabaseAdaption is the computation of tensor decomposition, and its complexity is𝒪​\(i​\(8​s3​d3​dc3​R\+6​s​d​dc​R2\)\)\\mathcal\{O\}\(i\(8s^\{3\}d^\{3\}d\_\{c\}^\{3\}R\+6sdd\_\{c\}R^\{2\}\)\)\. The calculation of FuturePrediction is𝒪​\(k3​lc\+k2​dc​ls\+k2​ls\+k​ls​d\)\\mathcal\{O\}\(k^\{3\}l\_\{c\}\+k^\{2\}d\_\{c\}l\_\{s\}\+k^\{2\}l\_\{s\}\+kl\_\{s\}d\)\. Thus the total computation is𝒪​\(i​\(8​s3​d3​dc3​R\+6​s​d​dc​R2\)\+k3​lc\+k2​dc​ls\+k2​ls\+k​ls​d\)\\mathcal\{O\}\(i\(8s^\{3\}d^\{3\}d\_\{c\}^\{3\}R\+6sdd\_\{c\}R^\{2\}\)\+k^\{3\}l\_\{c\}\+k^\{2\}d\_\{c\}l\_\{s\}\+k^\{2\}l\_\{s\}\+kl\_\{s\}d\)\.

## Appendix DLimitations and Future Work

One limitation of our approach is its reliance on the tensor decomposition via ALS \(Alternating Least Squares\), which is known to be sensitive to initialization and may converge to local optima\. However, in our streaming setting, the decomposition is only triggered when a significant mismatch between the current system and the new data is detected\. Since the data typically evolves smoothly, we initialize ALS using the previous decomposition result, which significantly reduces the number of iterations and helps maintain stability in practice\. Empirically, we observe that this warm\-start strategy leads to consistent convergence with minimal sensitivity\.

Another limitation is that our framework currently targets a mixture of linear time\-delay systems\. Extending it to nonlinear or hybrid systems is a promising direction for future work\. One possible approach is to integrate kernelized Markov parameters or to replace linear system approximations with local surrogate models, which can better capture nonlinearity in physical processes\.

Finally, although our current experiments focus on low\- to medium\-dimensional systems, scalability to high\-dimensional settings \(e\.g\., robotics with vision input or complex multi\-agent control\) remains an open challenge\. Investigating low\-rank tensor compression or structure\-aware decomposition techniques may help address this in future extensions\.

## Appendix EAdditional Results

![Refer to caption](https://arxiv.org/html/2605.26191v1/x5.png)Figure 5:Effectiveness of the results: Normalized spectral norms of Markov parameters estimated byDelayMixfor synthetic systems with input delays of 1 and 3 \(τ=1,3\\tau=1,3\)\. Each row corresponds to a Markov parameter \(i\.e\., a system response at a time step\), and each column within a block matrix represents a component of the parameter\. The magnitude of each entry reflects the influence of the input at that position\. Delays are revealed as leading low\-magnitude entries, allowingDelayMixto estimate the underlying system delays without prior knowledge\.Table 3:Prediction error \(mean±\\pmstd\) for each method on different datasets and lookback stepslsl\_\{s\}\.Boldis best,Underlineis second best\.This section provides additional empirical results that complement the experiments in Section[4](https://arxiv.org/html/2605.26191#S4)\. We first report aggregated forecasting performance across all benchmarks, then illustrate howDelayMixrecovers input–delay structures and tracks regime changes on synthetic data, and finally examine the sensitivity ofDelayMixto key hyperparameters\.

### E\.1Numerical results

Table[3](https://arxiv.org/html/2605.26191#A5.T3)reports the mean and standard deviation of forecasting errors across all datasets forDelayMixand its baselines, with the best and second\-best results shown inboldandunderline, respectively\. As discussed in Section[4](https://arxiv.org/html/2605.26191#S4), our method outperforms state\-of\-the\-art approaches\.

### E\.2Effectiveness ofDelayMix

Fig\.[5](https://arxiv.org/html/2605.26191#A5.F5)shows the normalized spectral norms of the Markov parameters estimated byDelayMixon synthetic datasets generated from two systems with input delaysτ=1\\tau=1andτ=3\\tau=3\. Each row corresponds to a time step in the system’s response, and each column in the block matrix represents a Markov parameter component\. The spectral norm of each column reflects the strength of the input’s influence on that component\.

A system delay appears as a sequence of zeros or low\-magnitude values at the beginning \(left side\) of the matrix\. As shown,DelayMixaccurately identifies delays of 1 and 3 by detecting these characteristic patterns without any prior knowledge\. This illustrates its ability to recover input delay structures and adaptively model a mixture of time\-delay systems through estimated Markov parameters\. While the delay patterns are clearly visible in this experiment, their clarity may be reduced under high noise, overlapping delays, or frequent switching\. Nevertheless, we observe that they remain identifiable in practice\.

### E\.3Sensitivity to Hyperparameters

We evaluated the sensitivity ofDelayMixto two key hyperparameters \(RRandρ\\rho\)on the Ship\-OOD dataset with forecasting horizon fixed tols=30l\_\{s\}=30\. Figure[6](https://arxiv.org/html/2605.26191#A5.F6)reports the mean MAE and MSE when varyingR∈\{2,4,8,10\}R\\in\\\{2,4,8,10\\\}andρ∈\{0\.5,0\.6,0\.7,0\.8,0\.9,1\.0\}\\rho\\in\\\{0\.5,0\.6,0\.7,0\.8,0\.9,1\.0\\\}\. Both metrics change only mildly across the explored ranges, and no setting leads to a substantial degradation in performance, indicating thatDelayMixis robust to reasonable choices ofRRandρ\\rho\. We observe slightly better errors for intermediate values, but the overall trends remain flat, suggesting thatDelayMixdoes not require fine\-tuning of these hyperparameters to achieve strong forecasting performance on Ship\-OOD\.

![Refer to caption](https://arxiv.org/html/2605.26191v1/x6.png)Figure 6:Hyperparameter studies\.

## Appendix FExperimental Setting

Computing Infrastructure\.The configuration includes 2 \* Xeon Gold 6258R 2\.7Hz CPU, 12 \* 64GB DDR4 RAM \(768GB\), and NVIDIA RTX A6000 48GB GPU, which is sufficient for all the baselines\.

Metrics\.All models are typically assessed using cumulative mean\-squared errors \(MSE\) and mean\-absolute errors \(MAE\), which means the model’s performance is judged based on its total errors accumulated throughout the entire learning sequence\.

Hyperparameter Setting ofDelayMix\.We searched for the best results inρ∈\{0\.5,0\.6,0\.7,0\.8,0\.9,1\.0\}\\rho\\in\\\{0\.5,0\.6,0\.7,0\.8,0\.9,1\.0\\\}andR∈\{2,4,8,10\}R\\in\\\{2,4,8,10\\\}to minimize the errors in the validation dataset and fixedssto33\.

Baselines\.The baseline learning rate was selected from the range\[1​e−3,3​e−3,1​e−2,3​e−2\]\[1e\-3,3e\-3,1e\-2,3e\-2\]to minimize the errors in the validation data\. The other parameters are the same as those that we used the default parameters provided in the authors’ codes\.