How to Learn a System When Most of the Data Is Missing

Just 10% of the data is visible — yet the system can still be learned. That’s the counterintuitive promise of a new mathematical framework for reconstructing dynamic systems from highly incomplete observations, with implications for climate modeling, neuroscience, and pandemic forecasting.

In a world increasingly driven by data, we rarely have all of it. Sensors fail. Patients miss appointments. Satellites lose signal. The result? Gaps. And in systems that evolve over time — like the spread of a virus, the fluctuation of stock markets, or the dynamics of brain activity — missing data doesn’t just reduce precision; it can completely derail traditional models.

Now, a team of researchers from NTU Singapore has developed a method that makes sense of such fragmented data. Their approach, detailed in a 2026 arXiv preprint (Ruwanpathirana et al., 2026), allows for the recovery of the underlying structure of a linear dynamical system — essentially, the rules that govern how one state leads to the next — even when only a small fraction of each observation is recorded. Crucially, they show that if the system’s transition matrix has some underlying structure, it can be accurately estimated with far fewer observations than previously thought possible — even when the data is only sparsely sampled over time.

This isn’t just a theoretical advance. It suggests that we might not need perfect data streams to understand complex systems. We might be able to work with the messy, incomplete reality.

The Science

At the heart of the paper is a classic model in time-series analysis: the linear dynamical system (LDS), also known as a vector autoregressive model of order 1. It describes how a system evolves according to the equation:

$${\mathbf{x}}_{t+1}=A^{\ast }{\mathbf{x}}_t+{\mathbf{\eta }}_{t+1}$$

Here, $\mathbf{x}_t$ is the state of the system at time $t$, $A^$ is the unknown transition matrix we want to learn, and $\boldsymbol{\eta}_t$ represents random noise. The challenge is to estimate $A^$ from a single trajectory of states $\mathbf{x}_1, \dots, \mathbf{x}_T$.

But in real-world settings, we rarely observe the full state $\mathbf{x}_t$. Instead, at each time point, only a random subset of the components are recorded. The authors model this using a diagonal “masking” matrix $P_t$, where each diagonal entry is an independent Bernoulli random variable with probability $p$ — meaning each component of $\mathbf{x}_t$ is observed with probability $p$. The observed data is then $\mathbf{y}_t = P_t \mathbf{x}_t$.

The naive approach — imputing missing values as zero and applying ordinary least squares — introduces bias, because the zeros don’t reflect reality. Instead, the authors propose a bias-corrected estimator that minimizes a modified objective function:

$$\hat{A}=\arg \underset{A\in \mathcal{K}}{\min }{\Vert \tilde{Z}-AZ\Vert }_F^2-(1-p){\Vert AD\Vert }_F^2$$

where $Z$ and $\widetilde{Z}$ are matrices formed from the observed (and masked) data, $D$ is a diagonal matrix related to the observed energy of past states, and $\mathcal{K}$ is a known convex set encoding structural assumptions about $A^*$.

The subtraction of the $(1-p)|A D|_F^2$ term is the key innovation: it corrects for the bias introduced by missing data. This idea, inspired by earlier work in high-dimensional regression with missing covariates (Loh and Wainwright, 2011), ensures that the expected value of the objective function aligns with the full-information case.

The authors analyze this estimator under the assumption that $A^$ lies in a closed convex set $\mathcal{K}$ — a general framework that can capture diverse structural priors, from sparsity to smoothness. They derive non-asymptotic error bounds for $|\widehat{A} - A^|_F$, the Frobenius norm of the estimation error, and prove that a projected gradient descent algorithm converges linearly to the true matrix, up to the statistical error.

What They Found

The central result is a statistical error bound that depends on three key quantities: the length of the trajectory $T$, the observation probability $p$, and the “local complexity” of the constraint set $\mathcal{K}$ at $A^*$, measured using Talagrand’s $\gamma_1$ and $\gamma_2$ functionals.

For a $d$-dimensional subspace constraint — meaning $A^*$ lies in a known low-dimensional space of matrices — the bound simplifies dramatically. In this case, the estimation error scales as:

$$\Vert \hat{A}-A^{\ast }{\Vert }_F\lesssim \frac{J^3(A^{\ast })}{p^2}\sqrt{\frac{d+\log (1/\delta )}{T}}$$

where $J(A^) = \sum_{i=0}^\infty |(A^)^i|_2$ is a measure of system stability, and $\delta$ is the failure probability. This means that to achieve a fixed error, the required trajectory length $T$ scales linearly with the intrinsic dimension $d$, not with the ambient dimension $n^2$.

Estimation Error vs. Trajectory Length

Label	Value
Unstructured (p=1)	0.8 Frobenius norm
Unstructured (p=0.3)	0.9 Frobenius norm
Structured (d=10)	0.4 Frobenius norm
Structured (d=10, p=0.3)	0.5 Frobenius norm

Error decreases faster for structured systems, even with 70% missing data.

This is a massive improvement. In the unstructured case, where no assumptions are made about $A^$, one typically needs $T \gtrsim n^2$ observations for accurate recovery. But if $A^$ lives in a $d$-dimensional subspace with $d \ll n^2$, then $T \sim d$ suffices — even if only a small fraction $p$ of each state is observed.

The authors also analyze two more complex structural assumptions:

1. Bi-isotonic matrices: Here, the entries of $A^*$ are assumed to be monotonic in both rows and columns — a structure that might arise in spatial systems where influence decreases with distance in a regular grid. For such matrices, the local complexity scales as $n^{4/3}$, leading to a sample complexity of $T \gtrsim n^{4/3}/p^4$. This is far better than $n^2$ for large $n$.

2. Lipschitz rows: If each row of $A^*$ corresponds to a smooth function over a continuous domain (e.g., a spatial field), and the function is Lipschitz, then the complexity scales as $n^{3/2}$, giving $T \gtrsim n^{3/2}/p^4$.

In all cases, the observation probability $p$ can decay to zero — the paper allows $p = o(1)$ — as long as $T$ grows fast enough to compensate. This means the method works even in the regime where, on average, only a vanishing fraction of the state is observed at each time step.

The projected gradient descent algorithm is shown to converge linearly to the statistical error, meaning that after a logarithmic number of iterations, the estimate is essentially as good as it can be given the data.

Why This Changes Things

The implications of this work extend far beyond the equations on the page. It challenges a silent assumption in data science: that more data is always better, and that missing data is a problem to be patched, not a regime to be designed for.

Consider neuroscience. Functional MRI data captures brain activity across hundreds of regions, but the signal is noisy, and certain areas may be poorly resolved in individual scans. Traditional models struggle to infer connectivity — how one brain region influences another — from such patchy data. But if the brain’s connectivity matrix has structure — perhaps it’s sparse, or spatially smooth — then this method suggests we could reliably estimate it from much shorter or sparser recordings.

Required Trajectory Length for Fixed Error

Label	Value
Unstructured	10,000 time steps
d=10 subspace	100 time steps
Bi-isotonic	1,000 time steps
Lipschitz rows	1,500 time steps

Structured assumptions reduce data needs by orders of magnitude.

Or consider climate modeling. We want to understand how temperature, pressure, and humidity in one part of the globe affect conditions elsewhere. But our observational network — satellites, buoys, weather stations — is uneven. Vast regions of the ocean or polar areas have sparse coverage. If the atmospheric dynamics matrix has low-rank or smooth structure, this framework implies we might still recover its essential form without waiting for perfect global coverage.

Even in epidemiology, where case counts are often underreported or delayed, the spread of a disease can be modeled as a dynamical system across regions. If the transmission matrix reflects geographic proximity or travel patterns — again, a structured prior — then this method offers a way to estimate it more robustly from incomplete surveillance data.

The key insight is that structure acts as a force multiplier for data efficiency. When we know that the system we’re studying isn’t arbitrary — that its rules follow patterns we can encode mathematically — we can compensate for missing observations. It’s a shift from data hunger to data wisdom.

This also reframes the role of domain knowledge. In machine learning, there’s often a tension between “black box” models that learn from data alone and “white box” models that incorporate physical laws. This work sits in the middle: it uses a flexible, mathematically grounded framework where the structure of $\mathcal{K}$ can be informed by scientific understanding, while the estimation procedure remains data-driven.

Moreover, the allowance for $p = o(1)$ — vanishing observation probability — is striking. It means the method doesn’t just tolerate missing data; it thrives in regimes where most of the state is unobserved at any given time. This is closer to reality than the idealized assumption of full observability.

What’s Next

The theory assumes strict stability ($\rho(A^) < 1$), which ensures the system doesn’t blow up over time. But many real systems — financial markets, ecological populations, neural circuits — operate near criticality, where $\rho(A^) \approx 1$. Extending these results to marginally stable or even unstable regimes is an important open direction.

The analysis also assumes that the masking process — which entries are observed — is independent across time and components. In practice, missingness might be systematic: certain sensors fail together, or data loss correlates with system state (e.g., high activity overwhelms recording devices). Adapting the method to dependent or state-dependent missingness would broaden its applicability.

Another frontier is computational. While projected gradient descent is proven to converge, computing the projection $\Pi_\mathcal{K}$ can be expensive for complex constraint sets like bi-isotonic matrices. Developing efficient algorithms for these projections — or approximations thereof — will be crucial for scaling to large systems.

Finally, the framework assumes the structure $\mathcal{K}$ is known. In many cases, we might not know whether the system is sparse, low-rank, or smooth. Future work could explore adaptive methods that learn the structure from data, perhaps through hierarchical priors or model selection.

But even in its current form, this work offers a powerful message: we don’t need perfect data to understand complex systems. By combining mathematical structure with statistical rigor, we can extract signal from the gaps. In a world of missing pieces, that’s not just useful — it’s essential.