RobustLDS: Learning linear dynamical systems under adversarial corruptions

1 INTRODUCTION

• Robust estimation with dependent observations is less understood than adversarially robust regression with i.i.d. clean data.The target setting contaminates a fraction of the states in a single LDS trajectory.Source: §1 · pp. 1–2

1.2 Related work

• Prior robust-regression guarantees chiefly address i.i.d. clean data, while extending them to dependent observations is less understood.The introduction also surveys work on heavy-tailed noise, sparse process noise, and corrupted or partially observed states.Source: §1 · p. 1, §1.2 · pp. 2–3

1.1 Our contributions

• The paper introduces continuous relaxations of least-trimmed squares with alternating minimization, plus two group-sparse outlier estimators.The group-sparse formulations use either an ℓ1/ℓ2 penalty or hard constraints; the penalized estimator has non-asymptotic error rates.Source: §1.1 · p. 2
• The theoretical results are for an oblivious adversary and strictly stable systems, and the stated rates are likely near-optimal only when estimating both the system matrix and outliers.The authors describe these as the first theoretical results for robust LDS estimation under adversarial contamination from a single trajectory.Source: §1.1 · p. 2

1.3 Notation

• The notation section defines matrix row, column, and entry indices, common matrix norms, and shorthand for constant-factor bounds.It also specifies notation for identity matrices, Frobenius balls and spheres, and the Kronecker product.Source: §1.3 · p. 3

2 ROBUST ESTIMATORS BASED ON LEAST-TRIMMED SQUARES

• The estimator jointly selects uncontaminated adjacent observations and fits an A-constrained least-squares model, paralleling least trimmed squares.The outlier indicators are unknown, so the joint formulation is combinatorial and includes products of adjacent indicators.Source: §2 · p. 3
• The formulation assumes the number of contaminated observations, s, is known in advance.The paper follows prior robust-estimation work in making this assumption.Source: §2 · p. 3
• The meta-algorithm alternates between updating binary outlier indicators from a relaxed problem and refitting the system matrix with those indicators fixed.The relaxed indicator solution is rounded before the system-matrix update.Source: §2 · p. 4
• The proposed estimators use different continuous relaxations of the binary outlier indicators, including a sphere relaxation.For fixed system matrix, the sphere relaxation is a trust-region variant solvable with the Moré-Sorensen algorithm.Source: §2 · p. 4

2. Semi-definite program (SDP) relaxation.

• The SDP relaxation derives from a Shor relaxation of the indicator outer product, and for fixed A can be solved to obtain Y and z.The implementation uses an equivalent Burer-Monteiro method; its iteration count lacks theoretical bounds but is small in practice.Source: §2 · p. 4
• The bi-convex relaxation introduces b_t to represent whether adjacent indicators are both one, replacing their product with linear constraints.Its objective is bi-convex, and for fixed A the resulting problem is a linear program in b and z.Source: §2 · p. 4
• The convexified-Q relaxation shifts Q(A) to a positive-semidefinite matrix and relaxes z to the unit interval.Unlike the other relaxations described, it requires an estimate of A to be available.Source: §2 · pp. 4–5
• The relaxations are non-convex, and convergence is not known to hold for arbitrary random initialization.In the experiments, initialization uses A-constrained least squares on contaminated data.Source: §2 · p. 5

• To handle potentially unbounded outliers, the procedure filters adjacent state pairs and uses a diagonal matrix D to weight terms.Source: §3 · p. 5
• The authors use a simple filtering form; extending iterative filtering and weight learning to VAR is not straightforward because dependence breaks the underlying concentration arguments.Source: §3 · p. 5
• Two estimators jointly estimate A* and U* by exploiting the block-sparse structure of U*.Source: §3 · p. 6

1. Least squares with ℓ1/ℓ2 block-sparse penalization

• The first estimator solves a block-sparse penalized least-squares problem.Source: §1 · p. 6
• Because A is convex, the resulting program is convex.Source: §1 · p. 6
• The convex program can be solved efficiently in practice.Source: §1 · p. 6

2. Least squares with block-sparse hard constraints

• In the experiments, (P7) is solved using block coordinate descent and constrained weighted least squares via alternating minimization; (P8) uses alternating minimization.Source: §2 · p. 6
• The adversary chooses the outlier indices independently of the data, then may corrupt those observations arbitrarily.Source: §2 · p. 6
• If ρ(A*)<1, J(A*) is finite, though it can grow with n.Source: §2 · p. 6
• With probability at least 1−δ, any solution to (P7) satisfies the theorem’s error bound under Assumptions 1 and 2, for suitably large T.Source: §2 · p. 6
• The rate is described as near-optimal up to log factors for jointly estimating A* and U*, but suboptimal when estimating A* alone.Source: §2 · pp. 6–7
• For A-only estimation, the discussion says the s/T factor should be s²/T² and leaves adapting the proposed proof approach to dependent inliers as future work.Source: §2 · p. 7
• If s≤log(n/δ), T≳k² log⁴(n/δ) suffices; if s≥log(n/δ), T≳k²s log n log²(ns) suffices, with the quadratic dependence on k an artifact of the proof technique.Source: §2 · p. 7
• The proof sketch bounds two terms using concentration results, self-normalized martingale tails, and matrix-set arguments, then combines the bounds to obtain the theorem.Source: §2 · p. 7

4 EXPERIMENTS

• The experiments test random contamination by Cauchy outliers with location and scale 10, and report mean recovery error and variance across Monte Carlo simulations.For each fixed system, the authors run 20 simulations, repeat for 20 matrices, and initialize alternating minimization with A-constrained least squares on contaminated data.Source: §4 · pp. 7–8
• As T increases, error decreases for every estimator except Sphere-AM and OLS, while Sphere-AM can outperform Biconvex-AM when both ϵ and T are small.At ϵ=0, all estimators coincide with OLS; at ϵ=0.1, Biconvex-AM outperforms the other estimators, with LS-BSHC comparable.Source: §4 · p. 8
• At small ϵ, SDP-AM achieves smaller error, while Biconvex-AM performs best as ϵ increases; LS-BSHC remains comparable to Biconvex-AM across the tested range.OLS and Sphere-AM show the highest error overall, though Sphere-AM outperforms several alternatives when T=100 and ϵ≤0.02.Source: §4 · p. 8

AI use statement

• The authors state that they designed the algorithms and experiments, while generative AI assisted primarily with implementation code optimization and refactoring that they reviewed and tested.They report not using generative AI for the listed research and writing tasks, and accept responsibility for the final work.Source: AI use statement · p. 9

A TECHNICAL TOOLS

• The technical appendix recalls Talagrand’s γα functionals and majorizing-measure theorem, alongside concentration and tail-bound results used in the proof.The cited tools include results of Krahmer et al. (2014), Wainwright (2019), Abbasi-yadkori et al. (2011), and Laurent and Massart (2000).Source: A TECHNICAL TOOLS · pp. 13–15

B.1 Proof of Lemma 2

• Lemma 2 follows from a feasibility–optimality argument and decomposition over eS and eSc.The proof compares the optimal solution with the feasible pair (A*, U*).Source: B PROOF OF THEOREM 1 · p. 15, §B.1 · p. 21

B PROOF OF THEOREM 1

• The proof derives Theorem 1 by upper-bounding the right-hand side and lower-bounding the left-hand side of (B.3), then imposing a sufficient lower bound on T.The final bound follows by substituting the chosen r into (B.4) and simplifying.Source: B PROOF OF THEOREM 1 · pp. 15, 16, 19, 20

B.2 Proof of Lemma 3

• Lemma 3 applies Lemma 1 with an all-ones vector, takes a union bound over t = 2, . . . , T, and then takes a square root.Source: B PROOF OF THEOREM 1 · p. 16, §B.2 · p. 22

B.4 Proof of Lemma 5

• Lemma 5 uses the structure of uncontaminated data and standard concentration bounds to establish its stated high-probability control.The proof represents each state coordinate as a Gaussian linear form and applies a union bound over time and coordinates.Source: B PROOF OF THEOREM 1 · p. 16, §B.4 · pp. 23–24

B.5 Proof of Lemma 6

• Lemma 6 bounds a self-normalized martingale using sub-Gaussian innovations conditioned on the past, yielding the lemma’s high-probability statement.Source: B PROOF OF THEOREM 1 · pp. 16–17, §B.5 · p. 24

B.6 Proof of Lemma 7

• Lemma 7 establishes P(E3) ≥ 1 − δ by applying a matrix deviation theorem and using Gaussian-width results for eCK.The proof applies Theorem 2, identified as Theorem 3.1 of Krahmer et al. (2014).Source: B PROOF OF THEOREM 1 · p. 18

B.3 Proof of Lemma 4

• Lemma 4 bounds a Gaussian supremum using Lipschitz concentration and a union bound over coordinates, yielding P(E1a) ≥ 1 − δ.The argument relates the bound to Gaussian width.Source: §B.3 · pp. 22–23

B.6 Proof of Lemma 7

• The proof of Lemma 7 concludes that eSc ≥ T − 2s ≥ T/2, completing the proof.Source: §B.6 · p. 26

B.7 Proof of Lemma 8

• For Lemma 8, the proof applies Theorem 2 of Krahmer et al. (2014) and a union bound over t ∈ [T] to obtain a probability-at-least-1 − δ bound.Source: §B.7 · pp. 28–29

C.1 Robust regression under adversarial contamination

• Ordinary least squares is efficient on uncontaminated data, but adversarially replacing observations makes estimating β* a robust-regression problem.The review distinguishes inlier-outlier, Huber, and strong contamination models by how the outliers are selected and whether the adversary can observe the data.Source: §C.1 · p. 29

C.2 Robust learning of linear dynamical systems

• Fewer works address adversarial contamination in linear dynamical systems and higher-order VAR models than robust regression in the i.i.d. setting.Source: §C.2 · p. 31

C.2.1 Heavy-tailed covariates or process noise

• Heavy-tailed LDS and VAR studies use robust aggregation, regularized regression, or covariance-based estimators under assumptions such as multiple trajectories, stationarity, mixing, or finite moments.Examples include median aggregation across independent trajectories and robust or constrained Yule-Walker estimation.Source: §C.2.1 · pp. 31–32

C.3 Sparse process noise

• Kim and Lavaei (2026b) use coordinate-wise least absolute deviation followed by filtering and least squares, while Kim and Lavaei (2026a) analyze a Huber estimator for sparse adversarial process noise.For symmetric process noise, both works report an O(1/T) rate when the system is linear.Source: §C.3 · p. 32
• Least squares with an ℓ1 penalty on each innovative outlier recovers the system matrices under i.i.d. bounded-variance noise-plus-outlier and additional structural assumptions; Farahmand et al. (2011) study doubly robust dynamical-system smoothing.For that problem, they use a least-squares estimator with an ℓ1 penalty and give a provably convergent algorithm via coordinate descent.Source: §C.3 · p. 32

C.4 Learning from corrupted states

• Prior corrupted-state work includes multi-trajectory guarantees that do not apply to the single-trajectory setting and block-sparse thresholding for additive outliers under subset convexity and smoothness.Other cited settings include noisy observed states and entry-subsampled states.Source: §C.4 · p. 33

D ROBUST ESTIMATORS BASED ON LEAST-TRIMMED SQUARES

• Section D provides a detailed overview of the estimators introduced earlier.Source: D ROBUST ESTIMATORS BASED ON… · p. 33

D.1 Sphere relaxation

• The sphere relaxation begins with a binary selector z constrained to retain exactly T-s observations.Source: §D.1 · p. 33
• The sphere relaxation estimates z by solving a trust-region problem for fixed A, rounding to the T − s largest-magnitude entries, then estimating A by constrained weighted least squares.When z is fixed, the A step uses the A-constrained weighted least-squares problem.Source: §D.1 · pp. 33–34

D.2 Semi-definite program (SDP) relaxation

• The SDP relaxation replaces zz⊤ with a positive-semidefinite matrix Y; for fixed A it solves an SDP, then rounds z and estimates A by constrained weighted least squares.When Y = zz⊤ and z is binary, diag(Y) = z.Source: §D.2 · p. 34

D.3 Bi-convex relaxation

• The bi-convex relaxation introduces b_t = z_t z_(t+1), enforces the product through linear inequalities, and alternates convex estimation with rounding and constrained least squares.With A fixed, the relaxed program is a convex program with a linear objective; the rounded z values determine b before estimating A.Source: §D.3 · pp. 34–35

D.4 Convexifying Q(A) to estimate z

• Adding a shift α(A) = max(0, −λ_min(Q(A))) makes Q(A) positive semidefinite, enabling a convex relaxation for estimating z.The relaxed program replaces binary z with z in [0, 1]^T; the resulting estimate is rounded to retain the T − s largest-magnitude entries.Source: §D.4 · p. 35

E.1 Overview of the implementation

• Implementations use approximations when exact solvers are costly, including a tridiagonal trust-region solver for the sphere relaxation, Burer–Monteiro for the SDP, cvxpy for z steps, and FISTA for A steps.The sphere implementation uses SciPy’s solve_banded; the A-step FISTA projects onto the ℓ1 ball.Source: §E.1 · p. 36

• The experiments compare four contamination patterns: randomly placed Cauchy outliers, five-entry Cauchy bursts, randomly selected sign flips, and Gaussian noise on the largest-norm entries.The high-leverage condition contaminates the s entries with largest ℓ2 norm using Gaussian noise with standard deviation 8.Source: §E.2 · p. 37

E.2.1 Block burst contamination

• The plots vary T across 50, 100, 200, 400, and 800 at ϵ = 0, 0.02, and 0.1, and vary ϵ across five levels at T = 100, 400, and 800.Both plot sets use the block-burst contamination model.Source: §E.2.1 · pp. 37–38
• At ϵ = 0.1, ConvexQ(A)-AM has the lowest error, while Biconvex-AM and SDP-AM are comparable; as T increases, error falls for all except Sphere-AM and OLS.At small ϵ, SDP-AM initially leads, but Biconvex-AM reaches the lowest error as T grows; all estimators coincide with OLS at ϵ = 0.Source: §E.2.1 · p. 37, Figure 3 · p. 37
• Across T = 100, 400, and 800, OLS and Sphere-AM have the highest error; Biconvex-AM and ConvexQ(A)-AM perform best as ϵ increases.SDP-AM has lower error in the small-ϵ regime; Sphere-AM outperforms all estimators except SDP-AM when T ≤ 100 and ϵ ≤ 0.02.Source: §E.2.1 · p. 38, Figure 4 · p. 38

E.2.2 Sign flip contamination

• The plots vary T across 50, 100, 200, 400, and 800 at ϵ = 0, 0.02, and 0.1, and vary ϵ across five levels at T = 100, 400, and 800.Both plot sets use the sign-flip contamination model.Source: §E.2.2 · pp. 38–39
• Under sign flips, all estimators, including OLS, perform well; LS-BSP coincides with OLS, while Sphere-AM outperforms both and the remaining estimators show similar rates.Sign flips preserve entry magnitudes, and with the chosen λ, U = 0 is optimal, causing LS-BSP to collapse to OLS.Source: §E.2.2 · p. 38, Figure 5 · p. 38
• Across T = 100, 400, and 800, Biconvex-AM, SDP-AM, and LS-BSHC outperform the other estimators, while OLS and LS-BSP have the highest error and error rises with ϵ.Source: §E.2.2 · pp. 38–39, Figure 6 · p. 39

E.2.3 High-leverage contamination

• Across tested trajectory lengths, errors decay even with contamination, but OLS and LS-BSP decay more slowly; Biconvex-AM approaches A-constrained least squares on uncontaminated data.LS-BSP's choice of λ can absorb the added zero-mean Gaussian corruption, making it behave comparably to OLS in this setting.Source: §E.2.3 · p. 39, Figure 7 · p. 39
• LS-BSP's λ is chosen to dominate Gaussian innovations, so the added zero-mean Gaussian corruption can also be absorbed by λ.The authors attribute LS-BSP's OLS-like behavior in this setting to that choice.Source: §E.2.3 · p. 39
• As the contamination fraction increases, Frobenius error rises; OLS and LS-BSP have the highest error, while Biconvex-AM, SDP-AM, and LS-BSHC still achieve strong recovery.The reported pattern is similar across T = 100, 400, and 800.Source: §E.2.3 · pp. 39–40, Figure 8 · p. 40