Improved rates for prediction and identification for partially observed linear dynamical systems
Year: 2021

Abstract:
Identification of a linear time-invariant dynamical system from partial observa1 tions is a fundamental problem in control theory. Particularly challenging are sys2 tems exhibiting long-term memory. A natural question is how learn such systems 3 with non-asymptotic statistical rates depending on the inherent dimensionality (or4 der) d of the system, rather than on the possibly much larger memory length. We 5 propose an algorithm that given a single trajectory of length T with gaussian ob6 servation noise, learns the system with a near-optimal rate of ‹ O (︁» d T )︁ in H2 er7 ror, with only logarithmic, rather than polynomial dependence on memory length. 8 We also give bounds under process noise and improved bounds for learning a 9 realization of the system. Our algorithm is based on multi-scale low-rank approx10 imation: SVD applied to Hankel matrices of geometrically increasing sizes. Our 11 analysis relies on careful application of concentration bounds on the Fourier do12 main – we give sharper concentration bounds for sample covariance of correlated 13 inputs and forH∞ norm estimation, which may be of independent interest. 14

Identification of a linear time-invariant dynamical system from partial observa-1 tions is a fundamental problem in control theory. Particularly challenging are sys-2 tems exhibiting long-term memory. A natural question is how learn such systems3 with non-asymptotic statistical rates depending on the inherent dimensionality (or-4 der) 𝑑 of the system, rather than on the possibly much larger memory length. We5 propose an algorithm that given a single trajectory of length 𝑇 with gaussian ob-6 servation noise, learns the system with a near-optimal rate of ‹𝑂 (︁» 𝑑𝑇 )︁ in ℋ2 er-7 ror, with only logarithmic, rather than polynomial dependence on memory length.8 We also give bounds under process noise and improved bounds for learning a9 realization of the system. Our algorithm is based on multi-scale low-rank approx-10 imation: SVD applied to Hankel matrices of geometrically increasing sizes. Our11 analysis relies on careful application of concentration bounds on the Fourier do-12 main – we give sharper concentration bounds for sample covariance of correlated13 inputs and forℋ∞ norm estimation, which may be of independent interest.14
1 Introduction15
We consider the problem of prediction and identification of an unknown partially-observed linear16 time-invariant (LTI) dynamical system with stochastic noise,17
𝑥(𝑡) = 𝐴𝑥(𝑡− 1) + 𝐵𝑢(𝑡− 1) + 𝜉(𝑡) (1) 𝑦(𝑡) = 𝐶𝑥(𝑡) + 𝐷𝑢(𝑡) + 𝜂(𝑡), (2)
with a single trajectory of length 𝑇 , given access only to input and output data. Here, 𝑢(𝑡) ∈ R𝑑𝑢 are18 inputs, 𝑥(𝑡) ∈ R𝑑 are the hidden states, 𝑦(𝑡) ∈ R𝑑𝑦 are observations (or outputs), 𝜉(𝑡) ∼ 𝑁(0,Σ𝑥)19 and 𝜂(𝑡) ∼ 𝑁(0,Σ𝑦) are iid gaussian noise, and 𝐴 ∈ R𝑑×𝑑, 𝐵 ∈ R𝑑×𝑑𝑢 , 𝐶 ∈ R𝑑𝑦×𝑑, 𝐷 ∈ R𝑑𝑦×𝑑𝑢20 are matrices. Partial observability refers to the fact that we do not observe the state 𝑥(𝑡), but rather21 a noisy linear observation 𝑦(𝑡).22
As a simple and tractable family of dynamical systems, LTI systems are a central object of study for23 control theory and time series analysis. The problem of prediction and filtering for a known system24 dates back to [Kal60]. However, in many machine learning applications, the system is unknown and25 must be learned from input and output data. Identification of an unknown system is often a necessary26 first step for robust control [DMM+19, BMR18]. In a long line of recent work, the interplay between27 machine learning and control theory has borne fruit in an improved understanding of the statistical28 and online learning guarantees for prediction, identification, and control for unknown systems. In29 machine learning, LTI systems also serve as a simple model problem for learning from correlated30
Submitted to 35th Conference on Neural Information Processing Systems (NeurIPS 2021). Do not distribute.
data in stateful environments, and can give insight into understanding the successes of reinforcement31 learning [Rec19, TR19] and recurrent neural networks [HMR18].32
Partial observability poses a significant challenge to system identification: In the fully observed33 setting, given access to 𝑥(𝑡), there is no obstacle to learning the matrices directly through linear34 regression. However, in the partially observed setting, the most natural form of the optimization35 problem is non-convex.36
Systems exhibiting long-term memory are particularly challenging to learn. Restricting to marginally37 stable systems, this occurs when the spectral radius of 𝐴, 𝜌(𝐴), is close to 1, and it implies that the38 output at a particular time cannot be accurately estimated without taking into account inputs over39 many previous time-steps—on the order of 𝑂 Ä 1
1−𝜌(𝐴)
ä times steps. Such systems often arise in40
practice. A particular class of such systems are those exhibiting multiscale behavior, with different41 state variables that change on vastly different timescales [CR10]. For example, the body’s pH level is42 affected both by long-term changes on a timescale of days or weeks, as well as breathing rate which43 changes over a timescale of seconds. For such systems, it makes sense to discretize at the scale of the44 fastest changing variable, which leads to a long memory for the slowest-changing variable. With few45 exceptions, existing guarantees for learning partially observed LTI systems degrade as the memory46 length increases. However, counting the number of parameters in the model (1)–(2) suggests that the47 right measure of statistical complexity is the intrinsic dimensionality of the system, not the memory48 length. This leads to the following natural question.49
Question: How can we learn partially observed LTI systems with (non-asymptotic) statistical rates50 that depend on the intrinsic dimensionality of the system, rather than the memory length?51
Despite the simplicity of the question, little in the way of theoretical results are known. We focus52 on the particular problem of learning the impulse response (IR) function of the system—which fully53 determines its input-output behavior—in ℋ2 norm. This is a natural norm for prediction problems54 as it measures the expected prediction error under random input. Known guarantees for learning the55 IR depend on the memory length. One particularly undesirable consequence is that for a continuous56 system with time discretization ∆ going to 0, the memory scales as 1/∆ (while the system order57 stays constant), leading to suboptimal estimation by an arbitrarily large factor.58
Our key contribution is an algorithm and analysis that gives statistical rates that are optimal up to59 logarithmic factors. Unlike previous works, our rates depend on the system order 𝑑—the natural60 dimensionality of the problem—and only logarithmically on the memory length of the system. Our61 algorithm is based on taking a low-rank approximation (SVD) of the Hankel matrix, which is a62 widely used technique in system identification. We consider a multiscale version of this algorithm,63 where we repeat this process for a geometric sequence of sizes of the Hankel matrix. This is essential64 for obtaining a stronger theoretical guarantee. In the setting of zero process noise, we prove that our65
algorithm achieves near-optimal ‹𝑂 (︁»𝑑(𝑑𝑢+𝑑𝑦)𝑇 )︁ rates inℋ2 error for the learned system.66 Our analysis relies on careful application of concentration bounds on the Fourier domain to give67 sharper concentration bounds for sample covariance and ℋ∞ norm estimation, which may be of68 independent interest. While we consider our algorithm in a simple setting, we hope that this is a69 first step to understanding and improving more complex subspace identification algorithms. Indeed,70 SVD and related spectral methods are a standard step used in subspace identification algorithms71 such as N4SID; our analysis suggests that SVD has an important “de-noising effect”.72
We also give improved bounds for system identification, that is, learning the matrices 𝐴,𝐵,𝐶,𝐷73 using the Ho-Kalman algorithm [HK66], with ‹𝑂 (︁»𝐿𝑑(𝑑𝑢+𝑑𝑦)𝑇 )︁ rates.74 1.1 Notation75
Norms. We use ‖·‖ to denote the 2-norm of a vector. For a matrix 𝐴, let ‖𝐴‖ = ‖𝐴‖276 denote its operator norm, 𝜌(𝐴) denote its spectral radius (maximum absolute value of eigen-77
value), and ‖𝐴‖F denote its Frobenius norm. For a matrix-valued function 𝑀(𝑡) ∈ C𝑑1×𝑑2 ,78 ‖𝑀‖F := »∑︀ 𝑡 ‖𝑀(𝑡)‖ 2 F. Let 𝜎𝑟(𝐴) denote the 𝑟th singular value of 𝐴.79
Fourier transform. Given a matrix-valued function 𝐹 : Z → C𝑚×𝑛, define the (discrete-time)80 Fourier transform as the function “𝐹 : R/Z→ C𝑚×𝑛 given by “𝐹 (𝜔) = ∑︀∞𝑡=−∞ 𝐹 (𝑡)𝑒−2𝜋𝑖𝜔𝑡.81 Matrices. Given a sequence (𝐹 (𝑡))𝑎+𝑏−1𝑡=1 where each 𝐹 (𝑡) ∈ C𝑚×𝑛, define Hankel𝑎×𝑏(𝐹 ) as the82 𝑎𝑚 × 𝑏𝑛 block matrix such that the (𝑖, 𝑗)th block is [Hankel𝑎×𝑏(𝐹 )]𝑖𝑗 = 𝐹 (𝑖 + 𝑗 − 1). Given a83 sequence (𝐹 (𝑡))𝑎−1𝑡=0 where each 𝐹 (𝑡) ∈ C𝑚×𝑛, define the Toeplitz matrix as the block matrix such84 that the (𝑖, 𝑗)th block is [Toep𝑎×𝑏(𝐹 )]𝑖𝑗 = 𝐹 (𝑖 − 𝑗)1𝑖≥𝑗 . For a matrix 𝐴, let 𝐴⊤, 𝐴𝐻 , 𝐴† denote85 its transpose, Hermitian (conjugate transpose), and pseudoinverse, respectively. For a vector-valued86 function 𝑣 : {𝑎, . . . , 𝑏} → R𝑛, let 𝑣𝑎:𝑏 ∈ R(|𝑎−𝑏|+1)𝑛 denote the the vertical concatenation of87 𝑣(𝑎), . . . , 𝑣(𝑏).88
Control theory. For a matrix 𝐴 ∈ C𝑑×𝑑, define its resolvent as Φ𝐴(𝑧) = (𝑧𝐼 − 𝐴)−1. For a89 linear dynamical system 𝒟 given by (1)–(2), let Φ𝒟 = Φ𝑢→𝑦 denote the transfer function from90 𝑢 to 𝑦 (response to input). Then Φ𝒟 = Φ𝑢→𝑦 = 𝐶Φ𝐴𝐵 + 𝐷 = 𝐶(𝑧𝐼 − 𝐴)−1𝐵 + 𝐷. Let91 T := {𝑧 ∈ C : |𝑧| = 1} be the unit circle in the complex plane. For a matrix-valued function 𝐹 :92 T→ C𝑑1×𝑑2 , define theℋ2 andℋ∞ norms by93
‖Φ‖ℋ2 =   1
2𝜋 ∫︁ T ‖Φ(𝑧)‖2𝐹 𝑑𝑧 ‖Φ‖ℋ∞ = sup 𝑧∈T ‖Φ(𝑧)‖ .
For a function 𝐹 : N0 → C𝑑1×𝑑2 , define its Z-transform to be 𝒵[𝐹 ](𝑧) = ∑︀∞ 𝑛=0 𝐹 (𝑛)𝑧
−𝑛. Con-94 sidered as a function T→ C, we can take its ℋ2 and ℋ∞ norms. Overloading notation, we will let95 ‖𝐹‖ℋ𝑝 := ‖𝒵𝐹‖ℋ𝑝 for 𝑝 = 2,∞. Theℋ2 andℋ∞ norms can be interpreted as the Frobenius and96 operator norms of the linear operator from input to output, i.e., they measure the average power of97 the output signal under random or worst-case input, respectively. For background on control theory,98 see e.g., [ZDG+96].99
Constants. In proofs, 𝐶 may represent different constants from line to line.100
2 Main results101
We consider the problem of prediction and identification for an unknown linear dynamical sys-102 tem (1)–(2). Our main goal is to obtain error guarantees in ℋ2 norm, which determines prediction103 error under random input [OO19, Lemma 3.3].104 Problem 2.1. Consider the partially-observed LTI system 𝒟 (1)–(2) with gaussian inputs 𝑢(𝑡) ∼105 𝑁(0, 𝐼𝑑𝑢) for 0 ≤ 𝑡 < 𝑇 . Suppose that the system is stable, that is, 𝜌(𝐴) < 1, and that we observe106 a single trajectory of length 𝑇 started with 𝑥(0) = 0, that is, we observe 𝑢(𝑡) ∼ 𝑁(0, 𝐼𝑑𝑢) and 𝑦(𝑡)107 for 0 ≤ 𝑡 < 𝑇 .108
The goal is to learn a LTI system ‹𝒟 with the aim of minimizing ⃦⃦Φ‹𝒟 − Φ𝒟⃦⃦ℋ2 . Equivalently, letting109 𝐹 *(𝑡) = ® 𝐷, 𝑡 = 0 𝐶𝐴𝑡−1𝐵, 𝑡 ≥ 1 denote the impulse response function (also called the Markov parameters) of the system, the goal is110 to learn an impulse response ‹𝐹 minimizing ⃦⃦⃦𝐹 * − ‹𝐹 ⃦⃦⃦ ℋ2 = ⃦⃦⃦ 𝐹 * − ‹𝐹 ⃦⃦⃦ F .111
Note that learning 𝐹 * is sufficient to fully understand the input-output behavior of the system, but112 we may also ask to recover the system parameters 𝐴,𝐵,𝐶,𝐷 up to similarity transformation (see113 Theorem 2.3).114
Previous results [OO19, SRD19] roughly depend polynomially on the “memory” 11−𝜌(𝐴) , which115 blows up as the spectral norm of 𝐴 approaches 1. In the setting of zero process noise, our goal is116
to obtain rates that are ‹𝑂 Äpoly(𝑑,𝑑𝑢,𝑑𝑦)√ 𝑇 ä , with only poly-logarithmic dependence on 11−𝜌(𝐴) . See117 Figure 1 for a comparison.118
We assume that 𝜌(𝐴) < 1 because if 𝒟 is not stable, it is in general impossible to learn ‹𝒟 with119 finite ℋ2 error, as a system with infinite response can have arbitrarily small response on any finite120 time interval. However, it may still be possible to learn the response up to time 𝐿 ≪ 𝑇 in this121 case [SBR19]. The marginally stable case (𝜌(𝐴) = 1) is an important case we leave to future work.122
Method Rollout type Min # samples IR error
Least squares (IR) [TBPR17]
Multi 𝐿 𝜎 »
𝐿 𝑇 Least squares (IR) [OO19] Single 𝐿 𝜎 »
𝐿 𝑇 Nuclear norm minimization Multi min{𝑑2, 𝐿} 𝜎 »
𝐿 𝑇 [SOF20] Multi 𝑑 𝜎 »
𝑑𝐿 𝑇 rank-𝑑 SVD (Theorem 2.2) Single 𝐿 𝜎 »
𝑑 𝑇
Figure 1: Here, 𝐿 is the memory length for the system, which is ‹𝑂 Ä 11−𝜌(𝐴)ä for well-conditioned systems. Rollout type refers to whether we have access to a single trajectory or multiple trajectories. Min # samples refers to the minimum number of samples (up to log factors) before the bounds are operational. IR error refers to the error in the impulse response in Frobenius/ℋ2 norm. Logarithmic factors are omitted.
In our Algorithm 1, we first use linear regression to obtain a noisy estimate 𝐹 of the impulse re-123 sponse. Next, following standard system identification procedures, we form the Hankel matrix124 Hankel𝐿×𝐿(𝐹 ) with the entries of 𝐹 on its diagonals. Because the true Hankel matrix125
Hankel𝐿×𝐿(𝐹 *) =
à 𝐶𝐵 𝐶𝐴𝐵 · · · 𝐶𝐴𝐿−1𝐵
𝐶𝐴𝐵 𝐶𝐴2𝐵 ...
... . . . ... 𝐶𝐴𝐿−1𝐵 · · · · · · 𝐶𝐴2𝐿−1𝐵 í has rank 𝑑, we take a low-rank SVD 𝑅𝐿 of the Hankel matrix to “de-noise” the impulse response.126 We can then read off the estimated impulse response by averaging over the corresponding diagonal127 of 𝑅𝐿. For technical reasons, we need to repeat this process for a geometric sequence of sizes128 of the Hankel matrix: 𝐿 × 𝐿, 𝐿/2 × 𝐿/2, 𝐿/4 × 𝐿/4, and so forth. This is because the low-129 rank approximation objective for a ℓ × ℓ Hankel matrix encourages the diagonals that are Θ(ℓ) to130 be close—as those are the diagonals with the most entries—and hence estimates 𝐹 *(𝑡) well when131 𝑡 = Θ(ℓ). In other words, low-rank estimation for Hankelℓ×ℓ(𝐹 ) is only sensitive to the portion of132 the signal that is at timescale ℓ. Repeating this process ensures that we cover all scales.133
Our main theorem is the following.134
Theorem 2.2. There is a constant 𝐶1 such that following holds. In the setting of Problem 2.1,135 suppose that 𝐹 * is the impulse response function, 𝑇 is such that 𝑇 ≥ 𝐶1𝐿𝑑𝑢 log (︀ 𝐿𝑑𝑢 𝛿 )︀ , 𝜀trunc :=136 ⃦⃦
𝐹 *1[2𝐿,∞) ⃦⃦ ℋ∞ √ 𝑑𝑢 + ⃦⃦ 𝐺*1[2𝐿,∞) ⃦⃦ ℋ∞ ⃦⃦⃦ Σ 1/2 𝑥 ⃦⃦⃦ F , and 𝑀𝑥→𝑦 = (𝑂,𝐶,𝐶𝐴, . . . , 𝐶𝐴𝐿−1)⊤ ∈137 R(𝐿+1)𝑑×𝑑𝑦 . Let 0 < 𝛿 ≤ 12 and 𝜎 = » ‖Σ𝑦‖+ ‖Σ𝑥‖𝐿 log (︀ 𝐿𝑑𝑢 𝛿 )︀ ‖𝑀𝑥→𝑦‖2. Then with probabil-138
ity at least 1− 𝛿, Algorithm 1 learns an impulse response function ‹𝐹 such that139 ⃦⃦⃦‹𝐹 − 𝐹 *⃦⃦⃦ F = 𝑂 (︃ 𝜎   𝑑 (︀ 𝑑𝑦 + 𝑑𝑢 + log (︀ 𝐿 𝛿 )︀)︀ log𝐿 𝑇 + 𝜀trunc √ 𝑑 + ⃦⃦ 𝐹 *1(𝐿,∞) ⃦⃦ F )︃ .
In the absence of process noise (when Σ𝑥 = 𝑂), when 𝐿 and 𝑇 are chosen large enough, the first140 term dominates, and ignoring log factors, the dependence is 𝑂 (︁»
𝑑(𝑑𝑦+𝑑𝑢) 𝑇
)︁ . We expect this to be141
Algorithm 1 Learning impulse response through multi-scale low-rank Hankel SVD 1: Input: Length 𝐿 (power of 2), time 𝑇 . 2: Part 1: Linear regression to recover noisy impulse response 3: Let 𝑢(𝑡) ∼ 𝑁(0, 𝐼𝑑𝑢) for 0 ≤ 𝑡 < 𝑇 , and observe the outputs 𝑦(𝑡) ∈ R𝑑𝑦 , 0 ≤ 𝑡 < 𝑇 . 4: Solve the least squares problem
min 𝐹 :Supp(𝐹 )⊆[0,2𝐿−1] 𝑇−1∑︁ 𝑡=0 ‖𝑦(𝑡)− 𝐹 * 𝑢(𝑡)‖2 .
to obtain the noisy impulse response 𝐹 : [0, 2𝐿− 1] ∩ Z→ R𝑑𝑦×𝑑𝑢 . 5: Part 2: Low-rank Hankel SVD to de-noise impulse response 6: Let ‹𝐹 (0) = 𝐹 (0). 7: for 𝑘 = 0 to log2 𝐿 do 8: Let ℓ = 2𝑘. 9: Let 𝑅ℓ be the rank-𝑑 SVD of Hankelℓ×ℓ(𝐹 ) (i.e., argminrank(𝑅)≤𝑑 ‖𝑅−Hankelℓ×ℓ(𝐹 )‖).
10: For ℓ2 < 𝑡 ≤ ℓ, let ‹𝐹 (𝑡) be the 𝑑𝑦 × 𝑑𝑢 matrix given by ‹𝐹 (𝑡) = 1𝑡 ∑︀𝑖+𝑗=𝑡(𝑅ℓ)𝑖𝑗 , where (·)𝑖𝑗 denotes the (𝑖, 𝑗)th block of the matrix.
11: end for 12: Output: Estimate of impulse response ‹𝐹 . the optimal sample complexity up to logarithmic factors. However, in the presence of process noise,142 there is an undesirable factor of √ 𝐿 ‖𝑀𝑥→𝑦‖, which (for well-conditioned matrices) is expected to143
be 𝑂 Ä
1 1−𝜌(𝐴)
ä or 𝑂(𝐿). We leave it an open problem to improve the guarantees in this setting.144
Remark 1. The 𝐿-factor dependence on the process noise is unavoidable with the current algorithm:145 when the process noise has covariance Σ𝑦 = 𝐼 and decays after 𝐿 steps, it can cause perturbations of146 size 𝑂( √ 𝐿) compared to the noiseless system. Even in the case 𝑑 = 1, when the impulse response147
function is 𝑎𝑒−𝑘𝑡/𝐿 for a known 𝑘, the noise will cause the estimate of 𝑎 to be off by 𝑂( √ 𝐿), and148 hence the ℋ2 norm of the impulse response to be off by 𝑂(𝐿). Our algorithm only regresses on149 previous inputs, but in the presence of process noise, a better approach is to regress on both the150 previous inputs 𝑢(𝑡) and outputs 𝑦(𝑡) and then take a (weighted) SVD, as in N4SID [Qin06].151
Remark 2. A burn-in time of Ω(𝐿) is information-theoretically required to get poly(𝑑) rates. At-152 tempting to extrapolate an impulse response function from time 𝑜(𝐿) to time 𝐿 can magnify errors153 by exp(𝑑), because the finite impulse response of a system of order 𝑑 can approximate a polynomial154 of degree 𝑑− 1 on [0, 𝐿].155
We also show the following improved rates for learning the system matrices, by combining ℋ∞156 bounds for the learned impulse response with stability results for the Ho-Kalman algorithm [OO19].157 Because the input-output behavior is unchanged under a similarity transformation (𝐴,𝐵,𝐶) ← [158 (𝑊−1𝐴𝑊,𝑊−1𝐵,𝐶𝑊 ), we can only learn the parameters up to similarity transformation.159
Theorem 2.3. Keep the assumptions and notation of Theorem 2.2, suppose 𝒟 is observable and160 controllable, and let161
𝜀′ = 𝜎
  𝐿 (︀ 𝑑𝑦 + 𝑑𝑢 + log (︀ 𝐿 𝛿
)︀)︀ 𝑇 + 𝜀trunc.
Let 𝐻− = Hankel𝐿×(𝐿−1)(𝐹 *). Suppose that 𝜀′ = 𝑂(𝜎min(𝐻−)). Then with probability at least162 1−𝛿, the Ho-Kalman algorithm (Algorithm 1 in [OO19] with 𝑇1 = 𝐿, 𝑇2 = 𝐿−1) returns ̂︀𝐴, “𝐵, “𝐶163 such that there exists a unitary matrix 𝑊 satisfying164
max {︁⃦⃦⃦ 𝐶 − “𝐶𝑊 ⃦⃦⃦ F , ⃦⃦⃦ 𝐵 −𝑊−1“𝐵⃦⃦⃦ F }︁ = 𝑂(
√ 𝑑 · 𝜀′)⃦⃦⃦
𝐴−𝑊−1 ̂︀𝐴𝑊 ⃦⃦⃦ F = 𝑂
Ç 1
𝜎min(𝐻−) · √ 𝑑 · 𝜀′ · Ç ‖Φ𝒟‖ℋ∞ 𝜎min(𝐻−) + 1 åå .
As 𝐿 can be chosen to make 𝜀trunc negligible, this gives ‹𝑂 (︁»𝐿𝑑(𝑑𝑢+𝑑𝑦)𝑇 )︁ rates, however,165 with factors depending on the minimum eigenvalue of 𝐻 . This is an improvement over the166 ‹𝑂 (︁√𝑑 4»𝐿(𝑑𝑢+𝑑𝑦)𝑇 )︁ rates in [OO19].167 We prove Theorem 2.2 in Section 4 and Theorem 2.3 in Appendix B.168
3 Related work169
We survey two classes of methods for learning partially observable LDS’s, subspace identification170 and improper learning. With the exception of [RJR20], all guarantees have sample complexity171 depending on the memory length 𝐿, which we wish to avoid.172
3.1 Subspace identification173
The basic idea of subspace identification [Lju98, Qin06, VODM12] is to learn a certain structured174 matrix (such as a Hankel matrix), take a best rank-𝑘 approximation (using SVD or another linear di-175 mensionality reduction method), and learn the system matrices 𝐴,𝐵,𝐶,𝐷 up to similarity transfor-176 mation. Usage of spectral methods circumvents the fact that the most natural optimization problem177 for 𝐴,𝐵,𝐶,𝐷 is non-convex. However, classical guarantees for these methods are asymptotic.178
Recently, various authors have given non-asymptotic guarantees for system identification algo-179 rithms. [OO19] analyzed the Ho-Kalman algorithm [HK66] in this setting. [SRD19] consider the180 setting where system order is unknown and give an end-to-end result for prediction, while [TMP20]181 consider the problem of online filtering, that is, recovering 𝑥(𝑡)’s up to some linear transformation.182
An alternate, empirically successful approach is that of nuclear norm minimization or regulariza-183 tion [FPST13]. [SOF20] (building on [CQXY16]) give explicit rates of convergence, and show that184 the algorithm has a lower minimum sample complexity and is easier to tune.185
Our algorithm is based on the classical approach of taking a low-rank approximation of the Hankel186 matrix, but we repeat this process with Hankel matrices of sizes 𝐿×𝐿, 𝐿/2×𝐿/2, 𝐿/4×𝐿/4, and so187 forth; this is key modification that allows us to obtain better statistical rates. Our analysis builds on188 the analyses given in [OO19, SOF20]. As essential part of the analysis is analyzing linear regression189 for correlated inputs, where we extend the work of [DMR19] to MIMO systems, as explained below.190
3.1.1 Linear regression with correlated inputs191
An important step in obtaining non-asymptotic rates for system identification is analyzing linear192 regression for correlated inputs. The most challenging step is to lower-bound the sample covariance193 matrix of inputs to the linear regression. A lower bound, rather than a matrix concentration result, is194 sufficient [Men14, SMT+18, MT19]; however, a concentration result is obtainable in our setting.195
[TBPR17] give non-asymptotic bounds for learning the finite impulse response for a SISO system196 in ℓ∞ Fourier norm; however, they require 𝐿 rollouts of size 𝑂(𝐿) and hence Ω(𝐿2) timesteps.197 Addressing the more challenging single-rollout setting, [OO19] obtain bounds for a single rollout198 of ̃︀Ω(𝐿) timesteps, by using concentration bounds for random circulant matrices [KMR14] to de-199 rive concentration inequalities for the covariance matrix. These concentration inequalities for the200 covariance matrix were improved (by logarithmic factors) by [DMR19]. Although [DMR19] give201 an analysis in the SISO case, as we show in Theorem A.2, the results can be extended to the MIMO202 case with an 𝜀-net argument.203
3.2 Improper learning using autoregressive methods204
Instead of solving the statistical problem of identifying parameters, another line of work develops205 algorithms for regret minimization in online learning. The goal is simply to do well in predicting206 future observations, with small loss (regret) compared to the best predictor in hindsight; the learned207 predictor is allowed to be improper, that is, take a different functional form. In the stochastic case,208
this allows prediction almost as well as if the actual system parameters were known; however, the209 framework also allows for adversarial noise.210
One popular strategy for improperly learning the system is to learn a linear autoregressive filter over211 previous inputs and observations, or ARMA model. Naturally, because we are optimizing over a212 larger hypothesis class, the statistical rates depend on 𝐿 rather than the system order 𝑑.213
[GLS+20, Theorem 4.7] consider the problem of online prediction for a fully or partially observed214 LDS, and give a regret bound that depends polynomially on the memory length 𝐿. Their approach215 works even for marginally stable systems, that is, systems with 𝜌(𝐴) ≤ 1. See also [AHMS13,216 HSZ17, HLS+18, KMTM19, TP20, RJR20] for previous work using autoregressive methods.217
Of particular interest to us is [RJR20], which gives rates independent of spectral radius. Building218 on [HSZ17], they observe that it suffices to regress on previous inputs and outputs projected to a219 lower-dimensional space. Their algorithm works in the setting of process noise and competes with220 the Kalman filter, but only when 𝐴−𝐾𝐶 has real eigenvalues, where 𝐾 is the Kalman gain.221
4 Proof of main theorem222
In this section, we prove Theorem 2.2. The proof hinges on the following lemma, which shows223 that if we observe a low-rank matrix plus noise, then taking a low-rank SVD can have a de-noising224 effect, producing a matrix that is closer to the true matrix.225 Lemma 4.1 (De-noising effect of SVD). There exists a constant 𝐶 such that the following holds.226 Suppose that 𝐴 ∈ C𝑚×𝑛 is a rank-𝑟 matrix, ̂︀𝐴 = 𝐴 + 𝐸, and ̂︀𝐴𝑟 is the rank-𝑟 SVD of ̂︀𝐴. Then227 ⃦⃦⃦ ̂︀𝐴𝑟 −𝐴⃦⃦⃦
F ≤ 𝐶 √ 𝑟 ‖𝐸‖ . (3)
Compare this with the original error ⃦⃦⃦ ̂︀𝐴−𝐴⃦⃦⃦
F = ‖𝐸‖F, which can only be bounded by228 √︀
min{𝑚,𝑛} ‖𝐸‖. When applied to the 𝑑-SVD of the Hankel matrix, this gives a factor of √ 𝑑229
rather than √ 𝐿 for the error.230
Proof. We have231 ⃦⃦⃦ 𝐴𝑟 −𝐴 ⃦⃦⃦ F ≤ √ 2𝑟 ⃦⃦⃦ 𝐴𝑟 −𝐴 ⃦⃦⃦ 2
(4)
≤ √ 2𝑟 (︁⃦⃦⃦ 𝐴𝑟 −𝐴 ⃦⃦⃦ 2 + ⃦⃦⃦ 𝐴−𝐴 ⃦⃦⃦ 2 )︁ (5) ≤ 2 √ 2𝑟 ‖𝐸‖ (6)
where (4) follows from 𝐴𝑟 − 𝐴 having rank at most 2𝑟, (5) follows from the triangle inequality,232 and (6) follows from Weyl’s Theorem: ⃦⃦⃦ 𝐴𝑟 −𝐴 ⃦⃦⃦ 2 ≤ 𝜎𝑟+1(𝐴) ≤ 𝜎𝑟+1(𝐴) + ‖𝐸‖ = ‖𝐸‖.233
To prove Theorem 2.2, we will need to obtain bounds for 𝐹 : {0, 1, . . . , 2𝐿−1} → R𝑑𝑦×𝑑𝑢 learned234 from linear regression inℋ∞ norm. The following is our main technical result.235 Lemma 4.2. There are 𝐶1, 𝐶2 such that the following hold. Suppose 𝑦 = 𝐹 * * 𝑢 + 𝐺* *236 𝜉 + 𝜂 where 𝑢(𝑡) ∼ 𝑁(0, 𝐼𝑑𝑢), 𝜉(𝑡) ∼ 𝑁(0,Σ𝑥), 𝜂(𝑡) ∼ 𝑁(0,Σ𝑦) for 0 ≤ 𝑡 < 𝑇 ,237 and Supp(𝐹 *),Supp(𝐺*) ⊆ [0,∞). Let 𝐹 = argmin𝐹∈{0,...,𝐿}→R𝑑𝑦×𝑑𝑢 ∑︀𝑇−1 𝑡=0 |𝑦(𝑡) − (𝐹 *238
𝑢)(𝑡)|2, 𝑀𝐺* = (𝐺*(0), . . . , 𝐺*(𝐿))⊤ ∈ R(𝐿+1)𝑑×𝑑𝑦 , and 𝜀trunc = ⃦⃦ 𝐹 *1[𝐿+1,∞) ⃦⃦ ℋ∞ √ 𝑑𝑢 +239 ⃦⃦
𝐺*1[𝐿+1,∞) ⃦⃦ ℋ∞ ⃦⃦⃦ Σ 1/2 𝑥 ⃦⃦⃦ F . For 0 < 𝛿 ≤ 12 , 𝑇 ≥ 𝐶1𝐿𝑑𝑢 log (︀ 𝐿𝑑𝑢 𝛿 )︀ , −1 ≤ 𝑎 < 𝐿− 𝐿′,240 ⃦⃦
(𝐹 − 𝐹 *)1[𝑎+1,𝑎+𝐿′] ⃦⃦ ℋ∞
≤ 𝐶2
ñ… 1
𝑇
Ç  ‖Σ𝑦‖𝐿′ Å 𝑑𝑢 + 𝑑𝑦 + log Å 𝐿′
𝛿
ãã +   ‖Σ𝑥‖𝐿′𝐿𝑑𝑢 log Å 𝐿𝑑𝑢 𝛿 ã ‖𝑀𝐺*‖ å + 𝜀trunc ô with probability at least 1− 𝛿.241
In the case Σ𝑥 = 𝑂, this roughly says that the error in the learned impulse response, 𝐹 − 𝐹 *, over242 any interval of length 𝐿′, has all Fourier coefficients bounded in spectral norm by ‹𝑂 (︁»𝐿′(𝑑𝑢+𝑑𝑦)𝑇 )︁243 – what we expect if the error from linear regression is uniformly distributed over all frequencies.244
A complete proof is in Appendix A; we give a brief sketch. First, because the errors are Gaussian,245 the error from linear regression, 𝐹 − 𝐹 *, follows a Gaussian distribution. To bound its covariance,246 we lower-bound the smallest singular value of the sample covariance of the inputs (Lemma A.1,247 Appendix A.1). Here, the difficulty is that the inputs are correlated – the input at time 𝑡 is 𝑢𝑡:𝑡−𝐿.248 Fortunately, the translation structure means it is close to a submatrix of an infinite block Toeplitz249 matrix, which becomes block diagonal in the Fourier domain. This “decoupling” allows us to show250 concentration. Compared to the SISO setting in [DMR19], we require an extra 𝜀-net argument. Once251 we have a bound on the covariance, we can bound any ⃦⃦⃦ (◊ 𝐹 − 𝐹 *)(𝜔)⃦⃦⃦ by matrix concentration252
(Appendix A.2); to bound theℋ∞ norm it suffices to bound this over a grid of 𝜔’s (Lemma A.4).253
Bounding the error in ℋ∞ norm of the impulse response allows us to bound the error in operator254 norm of the Hankel matrix, as the following lemma shows.255
Lemma 4.3. For any 𝐹 : Z→ C𝑚×𝑛, we have ‖Hankel𝑎×𝑏(𝐹 )‖ ≤ ‖𝐹‖ℋ∞ .256
Proof. Note that when 𝑣 : Z → C𝑛, Supp(𝑣) ⊆ [0, 𝑏 − 1], we have (𝐹 * 𝑣)𝑏:𝑏+𝑎−1 =257 Hankel𝑎×𝑏(𝐹 )𝑣𝑏−1:0. Hence for any 𝑣 set to 0 outside of [0, 𝑏 − 1], using Parseval’s Theorem258 and the fact that the Fourier transform of a convolution is the product of the Fourier transforms, we259 have260
‖Hankel𝑎×𝑏(𝐹 )𝑣𝑏−1:0‖2 ≤ ‖𝐹 * 𝑣‖2 = ⃦⃦⃦“𝐹̂︀𝑣⃦⃦⃦
2 ≤ sup 𝜔∈[0,1] ⃦⃦⃦“𝐹 (𝜔)⃦⃦⃦ 2 ‖̂︀𝑣‖2 = ‖𝐹‖ℋ∞ ‖𝑣‖2
This shows that ‖Hankel𝑎×𝑏(𝐹 )‖ ≤ ‖𝐹‖ℋ∞ .261
Theorem 2.2 will follow from the following bound after an application of the triangle inequality.262 Lemma 4.4. There are 𝐶1, 𝐶2 such that the following holds for the setting of Problem 2.1.263 Suppose 𝐿 is a power of 2, and 𝑇 ≥ 𝐶1𝐿𝑑𝑢 log (︀ 𝐿𝑑𝑢 𝛿 )︀ . Let ⃦⃦ 𝐹 *1[𝐿+1,∞) ⃦⃦ ℋ∞ √ 𝑑𝑢 +264 ⃦⃦
𝐺*1[𝐿+1,∞) ⃦⃦ ℋ∞ ⃦⃦⃦ Σ 1/2 𝑥 ⃦⃦⃦ F and 𝑀𝑥→𝑦 = (𝑂,𝐶,𝐶𝐴, . . . , 𝐶𝐴𝐿−1)⊤ ∈ R(𝐿+1)𝑑×𝑑𝑦 . Then with265
probability at least 1− 𝛿, the output ‹𝐹 given by Algorithm 1 satisfies266 ⃦⃦⃦ (‹𝐹 − 𝐹 *)1[1,𝐿] ⃦⃦⃦
F ≤ 𝐶2
(︃  ‖Σ𝑦‖ 𝑑 (︀ 𝑑𝑦 + 𝑑𝑢 + log (︀ 𝐿 𝛿 )︀)︀ log𝐿
𝑇 +
  ‖Σ𝑥‖𝐿𝑑𝑑𝑢 log (︀ 𝐿𝑑𝑢 𝛿
)︀ 𝑇 ‖𝑀𝐺*‖+ 𝜀trunc √ 𝑑 )︃
Proof. We are in the situation of Lemma 4.2 with 𝐺*(𝑡) = 𝐶𝐴𝑡−11𝑡≥1. Let ℋℓ = Hankelℓ×ℓ(𝐹 )267 andℋ*ℓ = Hankelℓ×ℓ(𝐹 *). Suppose ℓ ≤ 𝐿 is even. Note that268
ℋℓ = Hankelℓ×ℓ(𝐹 *)⏟ ⏞ ℋ*ℓ + Hankelℓ×ℓ(𝐹 − 𝐹 *)
whereℋ*ℓ = Hankelℓ×ℓ(𝐹 *) is a rank-𝑑 matrix, with error term is bounded by269 ‖Hankelℓ×ℓ(𝐹 − 𝐹 *)‖ ≤ ⃦⃦ (𝐹 − 𝐹 *trunc)1[1,2ℓ−1] ⃦⃦ ℋ∞ by Lemma 4.3
< 𝐶
[︃… 1
𝑇
(︃  ‖Σ𝑦‖ ℓ Å 𝑑𝑢 + 𝑑𝑦 + log Å 𝐿′
𝛿 ãã +   ‖Σ𝑥‖ ℓ𝐿𝑑𝑢 log Å 𝐿𝑑𝑢 𝛿 ã ‖𝑀𝐺*‖+ 𝜀trunc ]︃ by Lemma 4.2 (7)
with probability at least 1− 𝛿. Let 𝑅ℓ be the rank-𝑑 SVD ofℋℓ. Then by Lemma 4.1,270 ‖𝑅ℓ −ℋ*ℓ‖F = 𝑂 Ä√ 𝑑 ‖Hankelℓ×ℓ(𝐹 − 𝐹 *)‖ ä . (8)
Now letting ‹𝐹 (𝑡) = 1𝑡 ∑︀𝑖+𝑗=𝑡(𝑅ℓ)𝑖𝑗 when ℓ2 < 𝑡 ≤ ℓ we have (using the fact that the mean271 minimizes the sum of squared errors)272
‖𝑅ℓ −ℋ*ℓ‖ 2 F ≥ ℓ∑︁ 𝑡= ℓ2+1 ∑︁ 𝑖+𝑗=𝑡 ‖(𝑅ℓ)𝑖𝑗 − 𝐹 *(𝑡)‖2F
≥ ℓ∑︁
𝑡= ℓ2+1
Å 𝑡 · ⃦⃦⃦‹𝐹 (𝑡)− 𝐹 *(𝑡)⃦⃦⃦2
F
ã ≥ Å ℓ
2 + 1 ã ℓ∑︁ 𝑡= ℓ2+1 Å⃦⃦⃦‹𝐹 (𝑡)− 𝐹 *(𝑡)⃦⃦⃦2 F ã .
Note that we only get a lower bound with a factor of ℓ if we restrict to 𝑡 that is Θ(ℓ), i.e., restrict to273 diagonals that have many entries. This is the reason we will have to repeat this process for multiple274 sizes. Hence275 ⃦⃦⃦
(‹𝐹 − 𝐹 *)1[ ℓ2+1,ℓ]⃦⃦⃦2F ≤ 1ℓ/2 ‖𝑅ℓ −ℋ*ℓ‖2F . Together with (8) and (7) this gives with probability ≥ 1− 𝛿 that276
⃦⃦⃦ (‹𝐹 − 𝐹 *)1[ ℓ
2 +1,ℓ] ⃦⃦⃦ F ≤ 𝐶
(︃  ‖Σ𝑦‖ 𝑑 (︀ 𝑑𝑦 + 𝑑𝑢 + log (︀ 𝐿 𝛿
)︀)︀ 𝑇 +
  ‖Σ𝑥‖𝐿𝑑𝑢 log (︀ 𝐿𝑑𝑑𝑢
𝛿 )︀ 𝑇 ‖𝑀𝐺*‖+ 𝜀trunc √ 𝑑√
ℓ
)︃
Replacing 𝛿 by 𝛿log2 𝐿 , using a union bound over powers of 2, and summing gives277 ⃦⃦⃦ (‹𝐹 − 𝐹 *)1[1,𝐿] ⃦⃦⃦
F = 𝑂
(︃  ‖Σ𝑦‖ 𝑑 (︀ 𝑑𝑦 + 𝑑𝑢 + log (︀ 𝐿 𝛿 )︀)︀ log𝐿
𝑇 +
  ‖Σ𝑥‖𝐿𝑑𝑑𝑢 log (︀ 𝐿𝑑𝑢 𝛿 )︀ log𝐿
𝑇 ‖𝑀𝐺*‖+ 𝜀trunc
√ 𝑑 )︃ .
278
Proof of Theorem 2.2. We have the bound in Lemma 4.4, and also the same bound279 for ⃦⃦⃦ (‹𝐹 − 𝐹 *)(0)⃦⃦⃦
F after applying Lemma 4.2 to (𝐹 − 𝐹 *)𝛿0. Finally, note that280 ⃦⃦⃦ (‹𝐹 − 𝐹 *)1(𝐿,∞)⃦⃦⃦ F = ⃦⃦ 𝐹 *1(𝐿,∞) ⃦⃦ F and use the triangle inequality.281
5 Experiments282
We compared three algorithms for learning the impulse response function: least-squares, and low-283 rank Hankel SVD with and without the multi-scale repetition. We include details of the experimental284 setup in Appendix D. Note that to reduce the number of scales, we consider use a slight modification285 of our Algorithm 1 which triples the size at each iteration instead.286
The plots show the error ⃦⃦ 𝐹 *1[1,𝐿] − 𝐹 ⃦⃦ 2 , where 𝐹 is the estimated impulse response on [1, 𝐿],287 averaged over 10 randomly generated LDS’s, as a function of the time 𝑇 elapsed. We consider288 systems of order 𝑑 = 1, 3, 5, 10, and memory lengths 𝐿 = 27, 81.289
Using SVD significantly reduces the error, supporting our theory which shows that SVD has a “de-290 noising” effect. Additionally, multiscale SVD has better performance than naive SVD when 𝑑 is291 moderate, 𝐿 is large, and data is limited, but the performance is similar in a data-rich setting.292
References293
[AHMS13] Oren Anava, Elad Hazan, Shie Mannor, and Ohad Shamir. Online learning for time294 series prediction. In COLT 2013 - The 26th Annual Conference on Learning Theory,295 June 12-14, 2013, Princeton University, NJ, USA, pages 172–184, 2013.296
[BMR18] Ross Boczar, Nikolai Matni, and Benjamin Recht. Finite-data performance guarantees297 for the output-feedback control of an unknown system. In 2018 IEEE Conference on298 Decision and Control (CDC), pages 2994–2999. IEEE, 2018.299
[BTR13] Badri Narayan Bhaskar, Gongguo Tang, and Benjamin Recht. Atomic norm denoising300 with applications to line spectral estimation. IEEE Transactions on Signal Processing,301 61(23):5987–5999, 2013.302
[CQXY16] Jian-Feng Cai, Xiaobo Qu, Weiyu Xu, and Gui-Bo Ye. Robust recovery of complex303 exponential signals from random gaussian projections via low rank hankel matrix re-304 construction. Applied and computational harmonic analysis, 41(2):470–490, 2016.305
[CR10] Shaunak Chatterjee and Stuart Russell. Why are dbns sparse? In Proceedings of306 the Thirteenth International Conference on Artificial Intelligence and Statistics, pages307 81–88. JMLR Workshop and Conference Proceedings, 2010.308
[DMM+19] Sarah Dean, Horia Mania, Nikolai Matni, Benjamin Recht, and Stephen Tu. On the309 sample complexity of the linear quadratic regulator. Foundations of Computational310 Mathematics, pages 1–47, 2019.311
[DMR19] Boualem Djehiche, Othmane Mazhar, and Cristian R Rojas. Finite impulse response312 models: A non-asymptotic analysis of the least squares estimator. arXiv preprint313 arXiv:1911.12794, 2019.314
[FPST13] Maryam Fazel, Ting Kei Pong, Defeng Sun, and Paul Tseng. Hankel matrix rank315 minimization with applications to system identification and realization. SIAM Journal316 on Matrix Analysis and Applications, 34(3):946–977, 2013.317
[GLS+20] Udaya Ghai, Holden Lee, Karan Singh, Cyril Zhang, and Yi Zhang. No-regret predic-318 tion in marginally stable systems. In COLT 2020 - The 33rd Annual Conference on319 Learning Theory, July 9-12, 2020, pages 1–44, 2020.320
[HK66] BL Ho and Rudolph E Kalman. Effective construction of linear state-variable models321 from input/output functions. at-Automatisierungstechnik, 14(1-12):545–548, 1966.322
[HLS+18] Elad Hazan, Holden Lee, Karan Singh, Cyril Zhang, and Yi Zhang. Spectral filtering323 for general linear dynamical systems. In Advances in Neural Information Processing324 Systems, pages 4634–4643, 2018.325
[HMR18] Moritz Hardt, Tengyu Ma, and Benjamin Recht. Gradient descent learns linear dynam-326 ical systems. The Journal of Machine Learning Research, 19(1):1025–1068, 2018.327
[HSZ17] Elad Hazan, Karan Singh, and Cyril Zhang. Learning linear dynamical systems via328 spectral filtering. In Advances in Neural Information Processing Systems, pages 1–2,329 2017.330
[Kal60] Rudolph Emil Kalman. A new approach to linear filtering and prediction problems.331 Journal of Basic Engineering, 82.1:35–45, 1960.332
[KMR14] Felix Krahmer, Shahar Mendelson, and Holger Rauhut. Suprema of chaos processes333 and the restricted isometry property. Communications on Pure and Applied Mathemat-334 ics, 67(11):1877–1904, 2014.335
[KMTM19] Mark Kozdoba, Jakub Marecek, Tigran T. Tchrakian, and Shie Mannor. On-line learn-336 ing of linear dynamical systems: Exponential forgetting in kalman filters. In Proceed-337 ings of the AAAI Conference on Artificial Intelligence, volume 33, pages 4098–4105,338 2019.339
[Lju98] Lennart Ljung. System identification: Theory for the User. Prentice Hall, Upper Saddle340 River, NJ, 2 edition, 1998.341
[Men14] Shahar Mendelson. Learning without concentration. In Conference on Learning The-342 ory, pages 25–39, 2014.343
[MT19] Nikolai Matni and Stephen Tu. A tutorial on concentration bounds for system identifi-344 cation. In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 3741–345 3749. IEEE, 2019.346
[OO19] Samet Oymak and Necmiye Ozay. Non-asymptotic identification of lti systems from347 a single trajectory. In 2019 American Control Conference (ACC), pages 5655–5661.348 IEEE, 2019.349
[Qin06] S Joe Qin. An overview of subspace identification. Computers & chemical engineer-350 ing, 30(10-12):1502–1513, 2006.351
[Rec19] Benjamin Recht. A tour of reinforcement learning: The view from continuous control.352 Annual Review of Control, Robotics, and Autonomous Systems, 2:253–279, 2019.353
[RJR20] Paria Rashidinejad, Jiantao Jiao, and Stuart Russell. Slip: Learning to predict in un-354 known dynamical systems with long-term memory. arXiv preprint arXiv:2010.05899,355 2020.356
[RV+13] Mark Rudelson, Roman Vershynin, et al. Hanson-wright inequality and sub-gaussian357 concentration. Electronic Communications in Probability, 18, 2013.358
[SBR19] Max Simchowitz, Ross Boczar, and Benjamin Recht. Learning linear dynamical sys-359 tems with semi-parametric least squares. arXiv preprint arXiv:1902.00768, 2019.360
[SGA20] Arnab Sarker, Joseph E Gaudio, and Anuradha M Annaswamy. Parameter estimation361 bounds based on the theory of spectral lines. arXiv preprint arXiv:2006.12687, 2020.362
[SMT+18] Max Simchowitz, Horia Mania, Stephen Tu, Michael I Jordan, and Benjamin Recht.363 Learning without mixing: Towards a sharp analysis of linear system identification. In364 Conference On Learning Theory, pages 439–473, 2018.365
[SOF20] Yue Sun, Samet Oymak, and Maryam Fazel. Finite sample system identification: Im-366 proved rates and the role of regularization. In 2nd L4DC 2020 - Learning for Dynamics367 & Control, June 11-12, 2020, 2020.368
[SRD19] Tuhin Sarkar, Alexander Rakhlin, and Munther A Dahleh. Finite-time system369 identification for partially observed lti systems of unknown order. arXiv preprint370 arXiv:1902.01848, 2019.371
[TBPR17] Stephen Tu, Ross Boczar, Andrew Packard, and Benjamin Recht. Non-asymptotic372 analysis of robust control from coarse-grained identification. arXiv preprint373 arXiv:1707.04791, 2017.374
[TMP20] Anastasios Tsiamis, Nikolai Matni, and George Pappas. Sample complexity of kalman375 filtering for unknown systems. In Learning for Dynamics and Control, pages 435–444.376 PMLR, 2020.377
[TP20] Anastasios Tsiamis and George Pappas. Online learning of the kalman filter with378 logarithmic regret. arXiv preprint arXiv:2002.05141, 2020.379
[TR19] Stephen Tu and Benjamin Recht. The gap between model-based and model-free meth-380 ods on the linear quadratic regulator: An asymptotic viewpoint. In Conference on381 Learning Theory, pages 3036–3083. PMLR, 2019.382
[Ver18] Roman Vershynin. High-dimensional probability: An introduction with applications383 in data science, volume 47. Cambridge University Press, 2018.384
[VODM12] Peter Van Overschee and BL De Moor. Subspace Identification for Linear Systems.385 Springer Science & Business Media, 2012.386
[WJ20] Andrew Wagenmaker and Kevin Jamieson. Active learning for identification of linear387 dynamical systems. In COLT 2020 - The 33rd Annual Conference on Learning Theory,388 July 9-12, 2020, 2020.389
[ZDG+96] Kemin Zhou, John Comstock Doyle, Keith Glover, et al. Robust and optimal control,390 volume 40. Prentice hall New Jersey, 1996.391
Checklist392
1. For all authors...393
(a) Do the main claims made in the abstract and introduction accurately reflect the paper’s394 contributions and scope? [Yes]395
(b) Did you describe the limitations of your work? [Yes] The main limitations are the396 sensitivity to process noise (Remark 1 after Theorem 2.2, elaborated in Appendix E),397 the restriction to Gaussian noise, and the assumption 𝜌(𝐴) < 1.398
(c) Did you discuss any potential negative societal impacts of your work? [N/A] The399 work is primarily a theory result.400
(d) Have you read the ethics review guidelines and ensured that your paper conforms to401 them? [Yes]402
2. If you are including theoretical results...403
(a) Did you state the full set of assumptions of all theoretical results? [Yes]404 (b) Did you include complete proofs of all theoretical results? [Yes]405
3. If you ran experiments...406
(a) Did you include the code, data, and instructions needed to reproduce the main experi-407 mental results (either in the supplemental material or as a URL)? [Yes] Code is in the408 supplement. Executing either the Julia file or notebook directly produces the graphs.409
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they410 were chosen)? [Yes] Appendix D411
(c) Did you report error bars (e.g., with respect to the random seed after running exper-412 iments multiple times)? [No] Graphs are averages over 10 runs; visual inspection of413 raw data shows that the variability is small.414
(d) Did you include the total amount of compute and the type of resources used (e.g., type415 of GPUs, internal cluster, or cloud provider)? [No] Code runs in 5 minutes on laptop.416
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...417
(a) If your work uses existing assets, did you cite the creators? [N/A]418 (b) Did you mention the license of the assets? [N/A]419 (c) Did you include any new assets either in the supplemental material or as a URL? [N/A]420
421
(d) Did you discuss whether and how consent was obtained from people whose data422 you’re using/curating? [N/A]423
(e) Did you discuss whether the data you are using/curating contains personally identifi-424 able information or offensive content? [N/A]425
5. If you used crowdsourcing or conducted research with human subjects...426
(a) Did you include the full text of instructions given to participants and screenshots, if427 applicable? [N/A]428
(b) Did you describe any potential participant risks, with links to Institutional Review429 Board (IRB) approvals, if applicable? [N/A]430
(c) Did you include the estimated hourly wage paid to participants and the total amount431 spent on participant compensation? [N/A]432

Review summaries:
This work studies the estimation of linear dynamical systems with unobserved state. The authors propose a novel algorithm which relies on obtaining low-rank approximations to the Hankel matrix. They provide several results on the performance of their algorithm, including finite time bounds on the error it obtains in estimating the system’s impulse response. As compared to previous work, the main novelty of this paper is eliminating a factor of the horizon in the case where there is no process noise.

The paper proposes an algorithm for learning LTI systems with zero process noise. Given an input trajectory generated from independent Gaussian input, the algorithm estimates the impulse response function, which in turn can be used to estimate system parameters up-to similarity transformation.  The key idea is the multi-scale smoothing of the impulse response: A first estimate is obtained using linear regression. Denoising of the response function estimate is then performed by low-rank approximation of Hankel matrices of different sizes, with the insight being that the size of the matrix determines which parts of the response functions are preserved under the low-rank approximation.    The proposed algorithm achieves an error rate that scales only logarithmically with the system memory. The efficacy of the algorithm is demonstrated on a synthetic problem.

This work considers the problem of learning linear dynamical systems with partial observability both in terms of learning the impulse response and the system parameters upto similarity transformation (using Ho-Kalman algorithm). The main technical contribution is Algorithm 1 and Theorem 2.2 which presents the error bounds for Algorithm 1. This shows that the impulse response can be learned upto frobenius error $\sqrt{d/T}$, without any dependence on the mixing time/memory length of the dynamics. The main technical tool used in Algorithm 1 is the least squares estimation of the Hankel matrix at various scales to cover dependences from inputs at various time instants and denoising it with SVD since the true Hankel matrix is constrained to have rank d.   The experiments illustrate the advantage of the algorithm compared to least squares regression.