Online Learning of the Kalman Filter with Logarithmic Regret

In this paper, we consider the problem of predicting observations generated online by an unknown, partially observed linear system, which is driven by stochastic noise. For such systems the optimal predictor in the mean square sense is the celebrated Kalman filter, which can be explicitly computed when the system model is known. When the system model is unknown, we have to learn how to predict observations online based on finite data, suffering possibly a non-zero regret with respect to the Kalman filter's prediction. We show that it is possible to achieve a regret of the order of $\mathrm{poly}\log(N)$ with high probability, where $N$ is the number of observations collected. Our work is the first to provide logarithmic regret guarantees for the widely used Kalman filter. This is achieved using an online least-squares algorithm, which exploits the approximately linear relation between future observations and past observations. The regret analysis is based on the stability properties of the Kalman filter, recent statistical tools for finite sample analysis of system identification, and classical results for the analysis of least-squares algorithms for time series. Our regret analysis can also be applied for state prediction of the hidden state, in the case of unknown noise statistics but known state-space basis. A fundamental technical contribution is that our bounds hold even for the class of non-explosive systems, which includes the class of marginally stable systems, which was an open problem for the case of online prediction under stochastic noise.


Introduction
The celebrated Kalman filter has been a fundamental approach for estimation and prediction of timeseries data, with diverse applications ranging from control systems and robotics [Bertsekas, 2017, Durrant-Whyte andBailey, 2006] to computer vision [Coskun et al., 2017] and economics [Harvey, 1990]. Given a known system model with known noise statistics, the Kalman filter predicts future observations of a partially observable dynamical process by filtering past observations. When the underlying process is linear and the noise is Gaussian, the Kalman filter is optimal in the sense that it minimizes the mean square prediction error. Since Kalman's seminal paper [Kalman, 1960], the stability and statistical properties of the Kalman filter have been well studied when the system model is known.
Learning to predict unknown partially observed systems is a significantly more challenging problem. Even in the case of linear systems, learning directly the model parameters of the system results in nonlinear, non-convex problems [Yu et al., 2018]. Adaptive filtering algorithms address the problem of making observation predictions when the system model or the noise statistics are unknown or changing [Ljung, 1978, Moore and Ledwich, 1979, Lai and Ying, 1991, Ding et al., 2006. These adaptive filtering approaches are usually based on variations of extended least squares. Despite the importance of adaptive filtering in applications such as GPS, the regret of online filtering algorithms has not been considered in this classical literature.
In this paper, we consider the problem of predicting observations generated by an unknown, partially observable linear dynamical system in state-space form. We assume that the system dynamics and observation map are corrupted by Gaussian noise. Our goal is to find an online prediction algorithm that has provable regret bounds with respect to the Kalman filter that has access to the full system model. Our technical contributions are: System theoretic regret: We define a notion of regret that has a natural, system theoretic interpretation. The prediction error of an online prediction algorithm is compared against the prediction error of the Kalman filter that has access to the exact model, which is allowed to be arbitrary. Previous regret definitions [Kozdoba et al., 2019] required the model to lie in a finite set.
Logarithmic regret for the Kalman filter: We present the first online prediction algorithm with provable logarithmic regret upper bounds for the classical Kalman filter. In fact, we prove that with high probability the regret of our algoritm is of the order ofÕ(1), whereÕ hides poly log N terms, where N is the number of observations collected. Our algorithm has polynomial time complexity, requires linear memory, and is based on subspace system identification techniques [Qin, 2006]. Instead of optimizing over the state-space parameters, which is a non-convex problem, we convexify the problem by establishing an approximate regression between the next observation and past observations. Our analysis is based on the stability properties of the Kalman filter, tools for self-normalized martingales and matrices, and additional results for persistency of excitation developed in this paper.
Logarithmic regret for non-explosive systems: Our regret guarantees hold for the class of non-explosive systems, which includes marginally stable linear systems as well as as systems where the state can grow at a polynomial rate. This settles an open question and conludes that online prediction performance does not depend on the system stability gap 1 of the system. Although it was recently shown that the stability gap does not affect system identification [Simchowitz et al., 2019], whether the stability gap affects online prediction under stochastic noise was an open problem.
Regret analysis for other predictors: Our approach directly carries over to various interesting online predictors. For example, our analysis can be directly extended to the case of f −step ahead prediction of observations. Another extension focuses on the regret of hidden state predictors when the state-space basis representation is known a priori. The latter situation arises, for example, when the state-space model is known but the noise statistics are unknown. All these predictors enjoy similar logarithmic regret bounds.
Gap between model-free LQR and Kalman filter: One of the implications of our bounds is that learning to predict observations like the Kalman filter is provably easier than solving the online Linear Quadratic Regulator (LQR) problem, which in general requires O( √ N ) regret. In fact, recent results suggest that in the LQR case, the regret is lower bounded by Ω( √ N ) [Simchowitz and Foster, 2020]. This might not be surprising due to the fact that, in the absence of exogenous inputs, we cannot inject exploratory signals into the system.

Related work
Recently, there have been very important results addressing the regret of the adaptive Linear Quadratic Regulator (LQR) problem , Faradonbeh et al., 2017, Ouyang et al., 2017, Abeille and Lazaric, 2018, Dean et al., 2018, Mania et al., 2019, Cohen et al., 2019. The best regret for LQR is sublinear and of the order ofÕ( √ N ), where N is the numbers of state samples collected; an in-depth survey can be found in Matni et al. [2019]. When the system model is known, then the Kalman filter is the dual of the Linear Quadratic Regulator, suggesting that this duality can be exploited in deriving the regret of the Kalman filter. However, when the system model is unknown, the Linear Quadratic Regular and the Kalman filter are not dual problems [Tsiamis et al., 2019]. As the state is fully observed in LQR, the system identification in adaptive LQR reduces to a simple least squares problem. In the adaptive Kalman filter, the state is partially observed resulting in non-convex system identification problems requiring us to consider a different approach.
A related but different problem focuses on online prediction algorithms for systems without internal states (such as ARMA -autoregressive moving average) [Anava et al., 2013]. Prediction of observations generated by state space models in the case of exogenous inputs and adversarial noise but with a bounded budget was studied in Hazan et al. [2018]. The work closest to ours is the very recent work of Kozdoba et al. [2019], where regret bounds with respect to the Kalman Filter was studied for the first time but in the restricted context of scalar and bounded observations. The regret is shown to be linear, where the linear term is small but nonzero.
Our online algorithm is inspired by subspace identification techniques [Bauer et al., 1999]. The technical approach is based on classical results for the analysis of the least-squares estimator for time series [Lai and Wei, 1982], high-dimensional statistics [Vershynin, 2018] as well as modern results for finite sample analysis of system identification in both the fully observed [Faradonbeh et al., 2018, Simchowitz et al., 2018, Sarkar and Rakhlin, 2018 and the partially observed case [Hardt et al., 2018, Oymak and Ozay, 2018, Simchowitz et al., 2019, Sarkar et al., 2019, Tsiamis and Pappas, 2019.
Paper organization. In Section 2 we provide some background on the classical Kalman Filter and formulate the regret problem considered in this paper. In Section 3 we introduce the online learning algorithm while the regret guarantees are presented in Section 4. We conclude with generalizations and discussion of future work in Sections 5, 6. Detailed proofs can be found in the Appendix. Notation. With 2 we denote the Euclidean norm for vectors and the spectral norm for matrices. The spectral radius of a matrix A is denoted by ρ (A). The smallest singular value of a matrix A is denoted by σ min (A). By A * we denote the transpose of A. Unless explicitly stated, when using the standard O(N ),Õ(N ), o(N ) notation we hide all other quantities, e.g. system constants, system dimensions, logarithms of failure probabilities. TheÕ(N ) notation hides (powers of) logarithmic terms of N . The poly(x) notation means a polynomial function of x.

Problem Formulation
The Kalman filter considers the problem of predicting observations generated by the following statespace system: where x k ∈ R n is the state, y k , ∈ R m are the observations (outputs), A ∈ R n×n is the system matrix and C ∈ R m×n is the observation matrix. The time series w k , v k represent the process and measurement noise respectively and are modeled as zero mean i.i.d. Gaussian variables, independent of each other, with covariances Q and R respectively. The initial state is zero mean Gaussian with covariance Σ 0 and independent of the noises. The following assumption holds throughout this paper.
Let F k σ(y 0 , . . . , y k ) be the filtration generated by the observations y 0 , . . . , y k . Given the observations up to time k, the optimal predictionŷ k+1 at time k + 1 in the minimum mean square error (mmse) sense is defined as: In the case of system (1), the optimal predictor admits a recursive expression, known as the Kalman filter :x where e k y k − Cx k is the innovation noise process. Matrix K ∈ R n×m is called the Kalman filter gain, and can be computed based on A, C, Q, R-see (7) in Subection 2.1. Although the Kalman filter gives the optimal mmse prediction, it requires the system matrices A, C and noise covariances Q, R to be known. In this paper, we seek online learning algorithms that can predict observations based only on past observation data, without any knowledge of system matrices of noise covariances. To quantify the online prediction performance, we define the regret of our online learning algorithm with respect to the Kalman filter (3) that has full knowledge of system model (1). Our goal is to achieve sublinear regret, as defined in the following problem statement.
Problem 1. Assume that A, C, Q, R in system model (1) are unknown. Consider a sequence y 0 , y 1 . . . of observations generated by system (1). Letỹ k ∈ F k−1 be the prediction of an online learning algorithm based on the history y k−1 , . . . , y 0 andŷ k be the Kalman filter prediction (3) that has full knowledge of model (1). Define the regret: Fix a failure probability δ > 0. Our goal is to find a learning algorithm such that with probability at least 1 − δ: where o(N ) does not depend on δ.
Our regret definition has a natural system theoretic interpretation since it is defined with respect to the Kalman filter. In Section 5, we discuss an alternative regret definition.
In the following subsection we provide some background on the Kalman filter and specify some standard assumptions, which guarantee that the Kalman filter is well-defined.

Kalman Filter Background
The Kalman filter enjoys two critical properties, namely closed-loop stability and innovation orthogonality, that are now reviewed. The following standard assumption holds throughout the paper and guarantees that the Kalman filter is well-defined.
Assumption 2. The system matrix pair (A, C) is observable, i.e. the observability matrix: has rank n for all k ≥ n. The pair (A, Q 1/2 ) is controllable, i.e. the controllability matrix has rank n for all k ≥ n, and R is strictly positive definite.
The following result shows that under Assumption 2, the closed loop matrix A − KC of the Kalman filter is stable.
Proposition 1 (Anderson and Moore 2005). Consider system (1) under Assumption 2. The Kalman filter gain in (3) is computed by: where P is the positive definite solution to Moreover, the closed-loop matrix A − KC is stable, i.e. it has spectral radius ρ(A − KC) < 1.
Proposition 1 implies that the Kalman filter reaches steady state exponentially fast, allowing us to assume the following.
Assumption 3. We assume that the initial state covariance is Σ 0 = P , where P is defined in (7).
If Σ 0 = P , then we have to consider time-varying gains K k in (3). The condition Σ 0 = P guarantees that the Kalman filter (3) has stabilized to its steady-state so that the gain K is constant. Since the Kalman filter converges exponentially fast to its steady-state Anderson and Moore [2005], this is a very mild assumption; it is also standard Knudsen [2001].
The next assumption makes sure that system (3) is minimal.
Assumption 4. The pair (A, K) is controllable.
If the pair (A, K) is not controllable, then we can find a similarity transformationx k = Sx k such that: But sincex k = 0 this implies thatx 2,k = 0 for all k ≥ 0. Hence we could removex 2,k and consider a reduced system representation with onlyx 1,k+1 . The following assumption is for notational simplicity. It assumes that the largest eigenvalue of A − KC is simple.
Assumption 5. For some M > 0 and all t ≥ 0, the closed-loop matrix satisfies If the largest eigenvalue has larger multiplicity then we can just consider ρ(A − KC) + ǫ in the above bound, for sufficiently small ǫ.
In addition to the previous stability properties , the other nice property of the Kalman Filter is that the innovation sequence e k = y k −ŷ k is orthogonal (uncorrelated) and, by Gaussianity, also i.i.d. By the law of large numbers, this implies that the ℓ 2 accumulative error N k=0 y k −ŷ k 2 2 will be of the order of O(N ) almost surely. Predicting the true observations exactly is impossible in the stochastic noise setting, even if we know the system model.
Note that both systems (1), (3) can generate the same observations y k , i.e. the noise parameterization is not unique Van Overschee and De Moor [2012]. Another source of ill-posedness is that the state space parameterization is non-unique. Any similarity transformation S −1 AS, CS, S −1 QS − * generates the same observations. In the following section, we will address these problems by considering an alternative system representation.

Online Prediction Algorithm
The main idea of our online prediction algorithm is based on a system representation that has been used in the subspace system identification Bauer et al. [1999]. Let p be an integer that represents how far we look into the past. We define the vector of past observations at time k: Define also the matrix of closed-loop responses By expanding the Kalman filter (3) p-steps into the past, the observation at time k can be rewritten as Instead of optimizing over system parameters A, C, K, which results in a non-convex optimization problem, we optimize over (the higher dimensional) G p , which makes the problem convex. From

Algorithm 1 Online Prediction Algorithm
Input: β, λ, T init such that T init > β log T init for k = 0, . . . , T init − 1 do Observe y k ,ỹ k = 0 end for for i = 1, 2, . . . do end for end for an online learning perspective, this technique is also known as improper learning. Using this lifting, we can learn a least squares estimateG k,p by regressing outputs y t to past outputs Z t,p for t ≤ k: where λ > 0 is a regularization parameter. Then, to predict the next observation, we can compute: We could also use the recursive update: as long as the past p is kept constant. Due to the stability properties of the Kalman filter (Section 2.1), if we consider p past observations, then the bias term in equation (10) is of the order of ρ(A − KC) p x x−p 2 . Notice that for non-explosive systems the statex k−p can grow polynomially fast in the worst case. Even if x k−p remains bounded, keeping the past p constant would lead to a non-vanishing bias error (linear regret). Thus, to make sure that the prediction error decreases, we need to gradually increase the past horizon p. For this reason, inspired by the "doubling trick" Cesa-Bianchi and Lugosi [2006], we divide the learning in epochs, where every epoch is twice longer than the previous one. During every epoch we keep the past horizon constant. Since ρ(A − KC) p is exponentially decreasing, it is sufficient to slowly increase the past as p = O(log T ), where T is the epoch duration.
The pseudo-code of our online prediction approach can be found in Algorithm 1. Each epoch lasts from time T i , . . . , 2T i − 1, where i = 1, . . . , is the epoch, T i = 2 i−1 T init , and T init is a design parameter (the length of the first epoch). During every epoch, we keep the past p i = β log(T i ) constant, where β is a design parameter. Initially, from time 0 to T init − 1, we have a warm-up phase where we gather enough observations to start predicting. To make sure that p i < T i , we tune T init accordingly. Within an epoch, the least squares based predictor (12) can be implemented in a recursive way, which requires polynomial complexity and at most O(log T i ) memory. In the beginning of an epoch, when p i is updated, we re-initialize the recursive predictor based on the whole past, which requires polynomial complexity and O(T i ) memory. Hence, in total, after N collected samples, the computational complexity is polynomial and the memory requirement is O(N ). In Section 6, we discuss ways to modify the initialization when changing epochs without using the whole past, which can reduce the memory to O(log N ).
An important property of Algorithm 1 is that no knowledge about the dynamics or even the state dimension n is required. Note that there is a tradeoff between the bias error and statistical efficiency. Increasing the past horizon by selecting larger β leads to smaller bias error, but increases the sample complexity of learning G p since we have more unknowns; it is also harder to achieve persistency of excitation, i.e. to have a large enough smallest singular value ofV k .

Regret Analysis
In this section, we prove that with high probability the prediction regret is not only sublinear, but also of the order of poly log N (orÕ(1)), where N is the number of observations collected so far. The challenge in the non-explosive regime is that the observations grow unbounded polynomially fast (Ω( √ N)). Meanwhile, recent work in finite sample analysis of system identification Oymak and Ozay shows that the model parameters can be learned at a slower rate (O(1/ √ N )). Therefore these system identification results cannot be directly applied to obtain regret bounds for our problem. Nonetheless, we show that our online Algorithm 1 mitigates the effect of unbounded observations. As a result, the logarithmic regret bound ofÕ(1) remains valid even as we approach instability.
We provide two results, one for non-explosive systems (ρ(A) ≤ 1) and one for stable systems (ρ(A) < 1). Before we present the regret results, let us introduce some standard notions. Let a(s) = s d − a d−1 s d−1 · · · − a 0 be the minimal polynomial of matrix A, i.e. the minimum degree polynomial such that a(A) = 0. Denote its degree by d. We define the ℓ 1 norm of its coefficients as a 1 d−1 i=0 |a i |; the ℓ 2 norm a 2 is defined in a similar way. Let κ be the dimension of the largest Jordan block of A that is a associated with an eigenvalue on the unit circle (i.e. ρ(A) = 1). Let κmax be the largest Jordan block among all eigenvalues. In general, κ ≤ κ max ≤ d ≤ n.
Theorem 1 (Regret for non-explosive systems). Consider system (3) with ρ(A) ≤ 1. Let y 0 , . . . , y N be sequence of system observations withŷ 0 , . . . ,ŷ N being the respective Kalman filter predictions. Let y 0 , . . . ,ỹ N be the predictions of Algorithm 1 with and fix a failure probability δ > 0. There exists a N 0 = poly (n, β, κ, log 1/δ), independent of N , such that with probability at least 1 − δ, if N > N 0 then: whereÕ(1) hidesÕ(N 0 ) andÕ(T 2κ init ) terms. Theorem 1 provides the first logarithmic regret upper bounds for the general problem of Kalman filter prediction. The burn-in time N 0 is related to persistency of excitation conditions, i.e. initially we need enough samples to guarantee that the smallest singular value of the Gram matrixV k increases linearly with k. Our bounds do not depend on the stability gap 1/(1 − ρ(A)) and they do not degrade as we approach instability. However, they suggest, via β, that the stability radius ρ(A − KC) of the Kalman filter closed-loop matrix affects the difficulty of learning.
Interestingly, our bounds show that the problem of learning to predict observations like the Kalman filter is provably easier than the online LQR, in the case of unknown model. The latter requires in general regret of the order of √ N [Simchowitz and Foster, 2020]. This is another reason why the problems are not dual in the unknown model case. This gap might be expected since in the case of Kalman filter without exogenous inputs, there is no exploratory signal.
The upper bound also depends on the quantities d κmax and a 1 , both of which can be exponential in the dimension of the system state n in the worst case. This can happen, for example, if κ = n, i.e. the system is an n−th order integrator. Dependence of learning performance on the coefficients of the characteristic or minimal polynomial has been found in related work [Hardt et al., 2018]. This dependence can be improved in some cases-see for example the phase polynomial in Hazan et al. [2018], where there are no repeated eigenvalues. In our case, this dependence could perhaps be improved by applying the techniques of Simchowitz et al. [2019]. However, it is an open question whether it is possible to avoid the exponential dependence on κ, κ max . It might be possible that systems with long chain structure, e.g. integrators, are indeed harder to learn. In system theory it is known that even in the known model case, such systems might be difficult to observe. In open-loop system identification [Simchowitz et al., 2019], such a dependence also appears. It might be an inherent limitation of the problem, since fundamental quantities of the system, for example matrix A i or the observability matrix O i scale with i κ .
Both of the above issues are avoided in the case of stable systems (ρ(A) < 1), where we have the following result.
Theorem 2 (Regret for stable systems). Consider system (3) with ρ(A) < 1. Let y 0 , . . . , y N be sequence of system observations withŷ 0 , . . . ,ŷ N being the respective Kalman filter predictions. Let y 0 , . . . ,ỹ N be the predictions of Algorithm 1 with β = Ω 1 log(1/ρ(A−KC)) . Fix a failure probability δ > 0. Then there exists a N 0 = poly (n, β, log 1/ρ(A), log 1/δ) such that with probability at least 1 − δ, if N > N 0 then: Notice that for stable systems we no longer have quantities that depend exponentially on the dimension n. The main bound (15) does not depend on the stability gap 1/(1 − ρ(A)) However, via N 0 , the guarantees depend logarithmically on the inverse radius log 1/ρ(A). This quantity is related to the time needed for a stable system to approach stationarity.
The proofs of Theorem 1 and Theorem 2 can be found in the Appendix. In the next subsection, we provide an overview of the regret analysis and explain why the quantities d κmax and a 2 appear in the bound in Theorem 1. We also explain what changes in the case of stable systems addressed by Theorem 2.

Regret analysis overview
Recall the definition of the innovation error e k = y k −ŷ k . For brevity, we also define the error e k ỹ k −ŷ k between the online prediction of Algorithm 1 and the Kalman Filter prediction. By adding and subtractingŷ k in the first term, we obtain It is sufficient to prove that the square loss ℓ 2 : is logarithmic in N . Because the innovations are i.i.d., we have a martingale structure for the second term since e k ∈ F k , whileẽ k ∈ F k−1 . The martingale term will in general be small and can be bounded in terms of the square loss L N . In particular, the quantity is a self-normalized martingale and can be analyzed based on the techniques of Abbasi-Yadkori et al.
[2011], Sarkar and Rakhlin [2018], which imply that with high probability. Hence, we will focus on bounding the square loss L N . For the remaining section, we will assume that we are within one epoch i so that the past horizon p = p i and T = 2 i−1 T init are kept constant. For brevity, we omit the subscript p from all variables and write G, Then, the error between our online prediction and the Kalman filter prediction can be written as: The regression term is due to the noise e k entering the system. The truncation bias is due to using only p past observations and not all of them.
The key ingredients to analyze the cumulative error L N are i) the stability properties of the closed-loop matrix A − KC; ii) self-normalization properties of predictor (12); and iii) persistency of excitation for the past observations with high probability. By persistency of excitation we mean that the least singular value of the Gram matrixV k is increasing as fast as O(k) with high probability.
Regression term. We can rewrite the regression term as a product of two separate terms: The first term, S k−1V −1/2 k−1 is again a self-normalized martingale and can be analyzed based on the techniques of , Sarkar and Rakhlin [2018], which imply that the term grows logarithmically with T . The martingale property again comes from the fact that the innovation process e k of the Kalman Filter is i.i.d.-see Section 2.1.
The second term,V . . , Z p . It could be bounded using the following lemma which is inspired by Lai and Wei [1982].
Then, the following inequality holds: The intuition is that Z k−1 Z * k−1 appears inV k−1 and, hence, it cancels out the effect of Z k−1 . Unfortunately, we cannot directly use the above inequality forV k is not explicitly contained inV k−1 . However, we can exploit the fact that Z k does not change too fast compared to the most recent past Z k−1 , . . . , Z k−n .
Lemma 2 (ARMA-like representation). Let y 0 , y 1 . . . be observations generated by system (1). Fix a past horizon p and let a be the minimal polynomial of A with degree d. Then, the past observations satisfy the following recursion where Intuitively, the unbounded components of Z k are captured by the recent history Z k−1 , . . . , Z k−d and the residual δ k is bounded. Replacing Z k with (18) we obtain by two Cauchy-Schwarz inequality applications: The termsV −1/2 k−1 Z k−d+i in the sum are now indeed normalized and can be bounded using Lemma 1. ForV −1/2 k−1 δ k we exploit a new persistency of excitation result.
Lemma 3 (Uniform Persistency of Excitation). Consider the conditions of Theorem 1. Select a failure probability δ > 0. Let T = 2 i−1 T init for some fixed epoch i with p = β log T the corresponding past horizon. There exists a N 0 = poly(n, β, κ, log 1/δ) such that if T ≥ N 0 , then with probability at least 1 − δ: uniformly for all T ≤ k ≤ 2T − 1.
The above persistency of excitation condition holds uniformly over all times k as long as k ≥ N 0 . This is why the burn-in time N 0 appears in Theorem 1; if k is very small, then matrix k j=p Z j Z * j is not even invertible. A similar persistency of excitation result was proved in Tsiamis and Pappas [2019] for a fixed time k. However, the result of Lemma 3 is more general since it holds for all k.
Regularization and Truncation terms. For the regularization term we follow the same steps as with the regression one. Since matrix A − KC is stable, the truncation term decreases exponentially fast with p. System quantity κ governs how fast the observations grow polynomially. Parameter β should be large enough cancel out this polynomial rate. This explains why κ affects the choice of β in (13).
Stable Systems. If ρ(A) < 1, then we can exploit the fact that Z k converges exponentially fast to a stationary distribution. Hence the termV −1/2 k−1 Z k will effectively be self-normalized, without the need to express Z k as a function of the past observations. In particular, for stable systems we prove a new and stronger persistency of excitation result. Denote: Then, we have the following.
Lemma 4 (Uniform Persistency of Excitation: Stable case). Consider the conditions of Theorem 2. Select a failure probability δ > 0. Let T = 2 i−1 T init for some fixed epoch i with p = β log T the corresponding past horizon. There exists a N 0 = poly(n, β, log 1/ρ(A), log 1/δ) such that if T ≥ N 0 , with probability at least 1 − δ: uniformly for all T ≤ k ≤ 2T − 1.
Hence, the termV −1/2 k−1 Z k can be bounded by: where now the normalized term Γ −1/2 Z,k Z k behaves like a standard isotropic Gaussian variable. The term log 1/ρ(A) is due to the fact that it takes O(log 1/ρ(A)) time for the state to approach the stationary distribution (mixing time).

Extensions
In this section, we discuss generalizations of Algorithm 1 and the regret analysis.
Alternative regret definition Denote the system responses by be the sequence of system responses. Then, a parameterization for online prediction could beȳ g k = g 1 y k−1 + · · · + g k y 0 , for all k ≥ 0. Let G ρ,L g : g t 2 ≤ Lρ t be the set of system responses which decay exponentially for some L and ρ < 1, which are larger than M, ρ(A − KC) in Assumption 5. This set can include for example, stable IIR filters or FIR filters. Then, an alternative regret definition would be: The above definition captures the one in Kozdoba et al. [2019], where the unknown system lies in a finite set.
Since the observations increase at most polynomially fast and due to the properties of the Kalman filter, we can show that the differenceR N − R N depends on logarithmic terms of N . Hence our definition (4), which does not require any model restriction is general. The details can be found in the Appendix.
f -steps ahead predictor An immediate generalization of Algorithm 1 is to consider the f −steps ahead predictor, where f is some future horizon. Instead of predicting only the next observation, we predict the sequence y k , . . . , y k+f −1 . Denote the future observations and noises by: Similar to (8), we can establish a regression: , and T f is a lower triangular block Toeplitz matrix generated by I, CK, . . . , CA f −2 K. The optimal Kalman filter predictor in this case iŝ Hence, the regret can be defined as in (4), with the lowercase y replaced with uppercase Y . The online predictor (12) can be adapted here:Ỹ whereG k,f,p is obtained similar to (11) by regressing future observations Y t to past observations Z t from time p up to k − f . The logarithmic regret guarantees ofÕ(1) also hold with the final bound depending polynomially on f and T f 2 .
State prediction If we have some knowledge about the state, e.g. the state space basis and the state space dimension n, then we can use the f −step ahead predictor to predict the hidden statê x k . Notice that the Kalman filter state predictionx k can be rewritten as: If we know O f and the future horizon is large enough f ≥ n we can compute the state prediction: whereỸ k is our f −step ahead prediction and † denotes the pseudo-inverse. In this case the regret: will enjoy the same logarithmic guarantees. Hence, our algorithm can be used to solve the adaptive Kalman filter problem posed in Mehra [1970], Anderson and Moore [2005], where the dynamics A, C are known but the noise statistics Q, R are unknown, with logarithmic regret. If we do not know A, C, then we could estimate the range space of O f by performing singular value decomposition onG k,f,p . However, there are infinite representations O f S, for any invertible S, all of which can explain the same observations. The definition (23) is ill-posed since x k andx k might be based on different transformations S. Finding an alternative regret definition is subject of future work.
Logarithmic memory It is possible to achieve the logarithmic regret guarantees with logarithmic O(log N ) memory, by modifying the initialization step in the beginning T of every epoch i in Algorithm 1. For stable systems, we could just resetG T −1 andV T −1 to zero and λI respectively. This might not work for non-explosive systems, since Z T,p can be polynomially large in T . In this case, based on the regret analysis, we could initializeG T −1 ,V T −1 with the recent history Z T −1,p , . . . , Z T −d,p and y T −1 , . . . , y T −d , whered is an upper bound for the degree of the minimal polynomial a. This gives us control overV −1/2 T −1 Z T,p -see Section 4, and requires only O(log T ) memory.

Conclusion and Future Work
In this paper, we provided the first logarithmic regret upper bounds for learning the classical Kalman filter of an unknown system with unknown stochastic noise. Our regret analysis holds for nonexplosive systems and our bounds do not degrade with the system stability gap.
Going forward, our paper opens up several research directions. An open question that is whether we can define an appropriate regret notion in the case of state prediction, when matrices A, C are unknown, and prove logarithmic bounds. Another interesting direction is to study how the learning performance is affected by system theoretic properties, such as the exponential quantity d κ in the case of systems with long chain structure, e.g. κ-order integrators. Analyzing the regret of other online algorithms, e.g. extended least squares, is also an open problem. Another challenging problem for both prediction and system identification is the case of explosive systems. Although in the fully observed case, this problem has been studied [Faradonbeh et al., 2018, Sarkar andRakhlin, 2018], it remains open in the case of partially observable systems. Finally, in this work we considered that the state is only driven by stochastic noise. A more general problem to study is when we also have exogenous inputs. One of the challenges is that it is harder to prove persistency of excitation in the case of closed-loop systems.

A Notation and organization of the Appendix
Structure. In the first Sections B, C we review results from system theory and statistics. These include the main tools with which we will prove Theorems 1, 2. In Section D, we provide PAC bounds and persistency of excitation (PE) results for a fixed time k (pointwise) and fixed past horizon p. In Section E, we generalize those PAC bounds and PE results from pointwise to uniform over all times k in one epoch. In Section F we prove Lemma 1. By combining the uniform bounds and Lemma 1, we prove in Section G that the square loss suffered within one epoch is logarithmic with respect the length of the epoch. Hence, we can now prove Theorem 1-see Section H. In Section I, we analyze the case of stable systems and prove Theorem 2. Finally, in Section J, we show how the alternative regret definition (22) is equivalent to ours (4) up to logarithmic terms. Section K includes some technical results about logarithmic inequalities, which are used to show that the burn-in time N 0 depends polynomially on the various system parameters.
Notation. Before we proceed with the regret analysis, let us introduce some notation. A summary can be found in Table 1 We will analyze the performance of Algorithm 1 based mainly on a fixed epoch i. Since the past horizon p is kept constant during an epoch, we will drop the index p from Z k,p , G p ,G k,p ,V k,p and write Z k , G,G k ,V k instead. Similar to the past outputs Z k , we also define the past noises: The batch past outputs, batch past noises, and batch past Kalman filter states are defined as: This notation will simplify the sums k t=p Z t Z * t =Z kZ * k , k t=p E t Z * t =Ē kZk etc. Recall the definition of the correlations between the current innovation and the past outputs S k k t=p e t Z * t and the regularized autocorrelations of past outputsV k λI +Z kZ * k . The innovation sequence e k is i.i.d. zero mean Gaussian. Its covariance has a closed-form expression and is defined as:R Ee k e * k = CP C * + R, where P is the solution to the Riccati equation (7). Next we define the Toeplitz matrix T k , for some k ≥ 1: A useful property of system (3) is that the past outputs can be written as: The covariance of T p E t is denoted by: We define the covariance of the state predictions: and the covariance of the past outputs: Finally, let A = SJS −1 be the Jordan form of A. With the big O notation we also hide parameters like C 2 , K 2 , R 2 , S 2 , S −1 2 etc.  1. If the system is asymptotically stable ρ(A) < 1, then A i

For both cases
Proof. Proof of first part. By Gelfand's formula Horn and Johnson [2012], for every ǫ > 0, there The proof for the other system quantities is similar.
Proof of second part. Assume that A is equal to a n × n Jordan block corresponding to λ = 1. The proof for the other cases is similar. Then we have that: By Lemma B.2, we obtain: where the second inequality is classical, see Exercise 0.0.5 in Vershynin [2018].
Hence, we have: Proof of third part. In the general case of pairs (λ j , κ j ) of eigenvalues and Jordan block sizes, similar to the previous proof: In the above proofs we used a standard result for the norm of (block) Toeplitz matrices. A proof can be found in Tsiamis and Pappas [2019].
Lemma B.2 (Toeplitz norm). Let M ∈ R m 1 n×m 2 n , for some integers n, m 1 , m 2 be an (upper) block triangular Toeplitz matrix: where M i ∈ R m 1 ×m 2 , i = 1, . . . , n. Then: The following result is standard, for a proof see Tsiamis and Pappas [2019]. (3), with Γ k Ex kx * k . The sequence Γ k is increasing in the positive semi-definite cone.

Lemma B.3 (monotonicity). Consider system
The following results is also standard, but we include a proof for completeness. (3), with Γ k Ex kx * k . Assume that the system is stable with ρ(A) < 1. Then, the sequence Γ k satisfies: Γ k = AΓ k−1 A * + KRK * and converges to Γ ∞ ≻ 0, the unique positive definite solution of the Lyapunov equation:

Lemma B.4 (Lyapunov difference equation). Consider system
Proof. Since A is stable Γ ∞ = ∞ k=0 A k KRK * (A * ) k is well defined and solves the Lyapunov equation. Since (A, K) is controllable Γ ∞ is strictly positive definite: where the controllability matrix KR 1/2 AKR 1/2 . . . A n−1 KR 1/2 has full rank. It is unique since the operator L(M ) = M − AM A * is invertible; it has eigenvalues bounded below by 1 − ρ 2 (A).
Notice that Γ 0 = 0 Γ ∞ and by induction, we can show that Γ k Γ ∞ . Since Γ k is also monotone, it converges to the unique Γ ∞ . Now form the difference It is sufficient to find a τ such that: Since the norm of grows as fast as A τ 2 = O(ρ(A) τ −κmax+1 τ κmax ), it is sufficient to pick: By Lemma K.1, the order of τ is Lemma B.5 (Convergence rate: Lyapunov difference equation). Consider system (3), with Γ k Ex kx * k . Assume that the system is stable with ρ(A) < 1. Proof.

B.3 Proof of Lemma 2
We first prove a result for just the observations y k .
Proposition B.1 (ARMA-like representation). Consider system (3) with a(s) = s d − a d−1 s d−1 − · · · − a 0 the minimal polynomial with degree d. Then the outputs can be represented by the ARMA recursion: where M s are matrices of appropriate dimensions such that: where a 1 denotes the ℓ 1 norm of the polynomial coefficients 1 + d−1 i=0 |a i |.
Proof. We start from the fact that: for t = 0, . . . , d. From the properties of the minimal polynomial we obtain that: The same will now hold for the past outputs: We can bound δ k by: Define: The fact that ∆ = O(d κmax−1 a 1 √ p) follows from the fact that A i−1 2 = O(i κmax−2 ).

C Statistical Toolbox
If k ≥ f 1 (p, δ) p + 128 mp 2 log 9 + p log 2 + p log 1 δ then with probability at least 1 − δ: Proof. Consider a 1/4-net N of the unit sphere S mp−1 . Then, from [Vershynin, 2018, Exercise 4.4.3 b]: Denote the partition of v in p blocks of length m as: Notice that U * v ∈ R k−p is zero-mean Gaussian with covariance matrix: where we define the convolutions: with r(0) = 1. Due to the Toeplitz structure: where the first inequality follows from Lemma B.2, while the second follows from the triangle inequality. By Cauchy-Schwartz: Moreover, for positive definite matrices the following bound holds Σ 2 F ≤ tr Σ Σ 2 = (k − p + 1)p. Hence, by Lemma C.2 below, we obtain Taking the union bound over the whole net N , we obtain where we used that |N | ≤ 9 mp Vershynin [2018]. From equation (C.2), we get that: k − p + 1 ≥ 128 mp 2 log 9 + p log 2 + p log 1 δ , which completes the proof.
The following result is standard but we include it to have a sense of the universal constants. It is a specialized version of the Hanson-Wright theorem Vershynin [2018].

C.2 Self-normalized martingales
The following theorem, which can be found in Tsiamis and Pappas [2019] is an extension of Theorem 1 in  and Proposition 8.2 in Sarkar and Rakhlin [2018].

Theorem C.1 (Cross terms,Tsiamis and Pappas [2019]). Let
for some integer r. Then, for any δ > 0, with probability at least 1 − δ, for all t ≥ 0 where M t,k ∈ R q×d . For some µ > 0 define: Fix a failure probability δ > 0. Then with probability at least 1 − δ: If EX k X * k is invertible for all k = 1, . . . , T , the result holds for µ = 0. Proof. Fix a k. By the following Lemma C.7 we obtain that we probability at least 1 − δ/T : The result follows by a simple union bound.
Lemma C.7. Consider v t ∈ R d ∼ N (0, I) i.i.d., for t = 1, . . . , k. Let X k ∈ R q be a linear combination: where M t ∈ R q×d . For someµ > 0 define: Then: If EX k X * k is invertible the result holds for µ = 0. Proof. An application of Jensen's inequality gives:

D PAC bounds and persistency of excitation for fixed-time and fixed-past
In this section, we include results for persistence of excitation and for identification of the system parameters for a fixed time instance k, fixed confidence δ and a fixed past horizon p. To avoid excess notation, we drop the dependence on p and k in this section.
Theorem D.1 (PAC bounds for identification). Consider system (3) with observations y 0 , . . . , y k . Fix a past horizon p and consider the notation of Table 1. Define k 1 (p, δ) p + 128 mp 2 log 9 + p log 2 + p log 1 With probability at least 1 − 5δ the following events hold at the same time: a) Persistency of excitation if k satisfies the following perstistency of excitation requirement Proof. Define the primary events: where matricesW t ,L t ,B t are appropriate Gram matrices that normalize the correlations: We will show that EZ and all of the above events occur with probability at least 1−δ each. Moreover Hence by a union bound: Part a: All primary events occur with probability at least 1 − δ. The fact that P (EX) ≥ 1−δ, P (EZ ) ≥ 1−δ follows by a Markov inequality argument-see Simchowitz et al. [2018]. The fact that P(E E ) ≥ 1 − δ follows from Lemma D.1. For the remaining events, we apply Theorem C.1. Notice thatR −1/2 e k and Σ −1/2 E E k are isotropic so that the conditions of Theorem C.1 hold.
Part b: Event E P E From Lemma D.2 below, we have that Part c: Event E cross We show that E cross ⊇ EZ ∩ E ZE . Conditioned on EZ ∩ E ZE , we have where we used the simplification δ mp(k−p+1) < 1 and log 5 m δ ≤ mp 2 log 3mp δ , for p ≥ 2.

D.1 Persistency of excitation proofs
First, we prove persistence of excitation of the past noises in finite time.
T p E k satisfy the conditions of Lemma C.1. Hence, under condition (D.10), with probability at least 1 − δ: Multiplying from both sides with Σ 1/2 Finally, from Tsiamis and Pappas [2019], we have Σ E σ R I.
Next, we prove persistency of excitation for the past outputs.
Lemma D.2 (Output PE). Consider the conditions of Theorem D.1 and the definition of E E , EX , E XE . If: The past outputs can be written as: For simplicity, we rewrite Σ −1/2 E T pĒk =Ū k , whereŪ k is defined similarly toĒ k but has unit variance components. As a result the sample-covariance matrix will be: The proof proceeds in two steps. First, we bound the cross-terms based on the events EX, E XE . Second, we show that if k − p + 1 is large enough, then the cross-terms are small enough so that Conditioned also on E XE : For simplicity denote: Let now u ∈ R mp , u 2 = 1 be an arbitrary unit vector. Then, consider the quadratic form

Conditioned on
we can bound the cross terms by:

Cross-terms are dominated
To complete the proof, it is sufficient to show that if (D.11) holds then also for any unit vector u. Define: To complete the proof, we apply the following Lemma D.3, where we exploit the fact that b ≥ σ min (Σ E ) ≥ σ R . It follows that it is sufficient to have: the following inequality is true: By optimizing over a, we obtain that the minimum over a is: It follows from the lemma below, which is more general.

F Proof of Lemma 1
We prove a slightly more general version.
Lemma F.1. Fix a p and consider the notation of Table 1. Consider an i ≥ 0. The following inequality is true: Proof. SinceV k is increasing in the positive semidefinite cone: Hence, it is sufficient to prove the inequality for i = 0. Recall thatV k =V k−1 + Z k Z * k . Consider the identity: where the last equality follows from the identity det(I + F B) = det(I + BF ). Rearranging the terms gives: where the inequality follows from the fact that the sequenceV k is increasing in the positive semidefinite cone and the elementary inequality: 1 − x ≤ log 1/x, for x ≤ 1.
Since the upper bound telescopes, we finally get

G Analysis within one epoch
We will analyze the ℓ 2 square loss for the duration of one epoch, from time T up to time 2T − 1 with fixed past horizon p = β log T . We have three cases: i) persistency of excitation is established T ≥ N 0 , where N 0 is defined in (E.3); ii) persistency of excitation is not established T < N 0 ; iii) warm-up epoch from 0 to T init − 1.
For now, we concentrate on the first two cases. Consider the ℓ 2 loss within the epoch: (G.1) Based on the notation of Table 1, the error between the Kalman filter prediction and our online algorithm is: Table 1. By Cauchy-Schwarz, the submultiplicative property of norm and by the fact that Z k−1V Hence bounding the L 2T −1 T consists of three steps, bounding the supremum , and the sum This is what we do in the following theorem. If moreover T ≥ N 0 then also: Proof. Uniform events occur with high probability. Consider the primary events E unif and E PE unif defined in (E.1), (E.2) and: In the following lemma, we bound term V −1/2 k−1 Z k 2 , which the key to obtaining bounds in the non-explosive regime. We will apply this in both cases i), ii).
. Consider the conditions of Theorem 1. Let a be the minimal polynomial of A with degree d and ∆ defined as in (B.4). Let T = 2 i−1 T init for some fixed epoch i with p = β log T the corresponding past horizon. Then, Proof. We replace Z k = a d−1 Z k−1 + · · · + a 0 Z k−d + δ k . Then by two applications of Cauchy-Schwarz: Now by Lemma 1 and Lemma 2 it follows that: The result follows from the fact that the sequenceV 2T −i is monotone.

G.2 Warm-up epoch
Proof. During this time, we haveỹ k = 0. Let Γ y,k = CΓ k C * +R. Then, By Lemma C.6 and by monotonicity of Γ k , with probability at least 1 − δ:

H Proof of Theorem 1
Recall that the regret can be decomposed in two terms: where L N is the square loss and the other term is a martingale. Square loss bound. Without loss of generality assume that N = 2T i − 1 = T init 2 i is the end of an epoch, where i is the total number of epochs. The number of epochs i depends logarithmically on N . Then the square loss L N can be written as: Then by Theorem G.1, Lemma G.2, and a union bound, with probability at least 1−(5 π 2 6 +1)δ −iδ 1 : L N =Õ(T 2κ init ) + poly(∆, a 2 2 , n, β, κ, log 1/δ, log 1/δ 1 ) Õ (N 0 ) +Õ (1) To complete the bound on L N , replace δ 1 with δ/i. Since i depends logarithmically on N we finally obtain that with probability at least 1 − (5 π 2 6 + 2)δ: L N =Õ(T 2κ init ) + poly(∆, a 2 2 , n, β, κ, log 1/δ) Õ (N 0 ) +Õ(1) Martingale term bound. Denote u k R −1/2 e k and z k R 1/2 (ŷ k −ỹ k ).
To apply Theorem C.1 we need to slightly modify the definition of the filtration. Let F t,i σ(F t ∪ {u t+1,1 , . . . , u t+1,i }), with F t+1 ≡ F t,m and define: By applying Theorem C.1 withF s we can bound the sum in terms of the square loss L N . With probability at least 1 − δ: where we used the fact that z *

I Stable case
In the case of stable systems, we can exploit the fact that the covariance matrix Γ k converges exponentially fast to a stead-state covariance Γ ∞ -see Lemma B.4. We have the following pointwise persistency of excitation result. Define the controllability matrix C t A t KR 1/2 . . . AKR 1/2 KR 1/2 (I.1) Since (A, K) is controllable we also have that rank(C t ) = n, for t ≥ n. As a result, the covariance matrix, is always strictly positive definite σ min (Γ t ) ≥ σ min (C n C * n ) > 0, for t ≥ n-see also Lemma B.4. Lemma I.1 (Stable: pointwise persistency of excitation). Consider system (3) with observations y 0 , . . . , y k . Pick a τ ≥ n. Define k 3 (τ, δ) τ + 128 mτ 2 log 9 + τ log 2 + τ log 1 With probability at least 1 − 3δ, if k ≥ k 3 (τ, δ), k 4 (k, τ, δ) The proof now is similar to Theorem D.1 and Lemma D.2. We will show that the cross terms are dominated. Define the primary events: where matricesW t , W are: The events EX, E XE occur with probability at least 1 − δ each-see proof of Theorem D.1 and Theorem C.1. By Lemma C.1, if k ≥ k 3 (τ, δ), then also E E occurs with probability at least 1 − δ. By a union bound EX , E E , E XE occur with probability at least 1 − 3δ. What remains to show is that these three events imply (I.3). The remaining proof is omitted since it is identical with the one of Lemma D.2.

I.1 Proof of Lemma 4
It follows from the lemma below, which is more general.

J Alternative regret definition
In this section, we sketch how the online learning definition (22) of regret, i.e. the best linear predictor if we knew all N data beforehand, is equivalent to our definition (4).
Lemma J.1. Consider system (3) with ρ(A) ≤ 1. Let y 0 , . . . , y N be sequence of system observations withŷ 0 , . . . ,ŷ N being the respective Kalman filter predictions. Fix a failure probability δ > 0. There exists a N 0 = poly (log 1/δ) such that with probability at least 1 − δ, if N > N 0 then: Z k = y * k−p · · · y * k−1 , for k = 0, . . . , N with y t = 0 if t < 0. Define E + e 0 · · · e N , Y + y 0 · · · y N , Z Z 0 · · · Z N , LetG = g 1 · · ·g p be the solution of min g 1 ,...,gp N k=0 y k − p t=1 g t y k−t 2 . We can show that:G = Y Z * (ZZ * ) −1 . Hence, where · F denotes the Frobenius norm. But we have: Replacing Y in (J.5), we obtain: But the term E 2 F = N k=0 y k −ŷ k 2 2 is the Kalman filter prediction error. Using Theorem D.1, we can show that there exists a N 0 such that ZZ * (k − p + 1)σ R I with high probability if N ≥ N 0 . Meanwhile, using Theorem C.1 combined with persistency of excitation, we can show that with probability at least 1 − δ: k ≥ c log k Then, a sufficient condition for the above inequality to hold is: Proof. If c ≤ e, then the inequality is satisfied for all k > 0. To see why this holds consider f (k) = k − e log k. The minimum is attained at f (e) = e − e log e = 0. Hence, k ≥ e log k ≥ c log k.
Next, we analyze the case c > e. We have that the function k − c log k is increasing for k ≥ c. Moreover, 2c log 2c ≥ c. As a result if k ≥ 2c log 2c then also: k − c log k ≥ 2c log 2c − c log(2c log 2c) = c log 2c − c log log 2c ≥ c log 2c − c e log 2c ≥ 0 where we used Lemma K.3.
Lemma K.2. Let c be a positive constant. Consider the inequality: k ≥ c log 2 k Then, a sufficient condition for the above inequality to hold is: k ≥ max 4c log 2 4c, 4c log 4c, 1 Proof. If c ≤ 1, then the inequality is satisfied for k ≥ 1. To see why this holds define f (k) = k − log 2 k. Its derivative f ′ (k) = 1 − 2 log k k is always positive for k ≥ 1 since from the proof of Lemma K.1 k ≥ 2 log k. Hence f (k) ≥ f (1) = 1.