extfieldsofstudy:
Computer Science
Engineering
Mathematics

s2fieldsofstudy:
Computer Science
Mathematics

text:
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.


Deep Stacked Ensemble Learning Model for COVID-19 Classification

: COVID-19 is a growing problem worldwide with a high mortality rate. As a result, the World Health Organization (WHO) declared it a pan-demic. In order to limit the spread of the disease, a fast and accurate diagnosis is required. A reverse transcript polymerase chain reaction (RT-PCR) test is often used to detect the disease. However, since this test is time-consuming, a chest computed tomography (CT) or plain chest X-ray (CXR) is sometimes indicated. The value of automated diagnosis is that it saves time and money by minimizing human effort. Three significant contributions are made by our research. Its initial purpose is to use the essential finetuning methodology to test the action and efficiency of a variety of vision models, ranging from Inception to Neural Architecture Search (NAS) networks. Second, by plotting class activation maps (CAMs) for individual networks and assessing classification efficiency with AUC-ROC curves, the behavior of these models is visually analyzed. Finally, stacked ensembles techniques were used to provide greater generalization by combining finetuned models with six ensemble neural networks. Using stacked ensembles, the generalization of the models improved. Furthermore, the ensemble model created by combining all of the finetuned networks obtained a state-of-the-art COVID-19 accuracy detection score of 99.17%. The precision and recall rates were 99.99% and 89.79%, respectively, highlighting the robustness of stacked ensembles. The proposed ensemble approach performed well in the classification of the COVID-19 lesions on CXR according to the experimental results.

Introduction
The coronavirus (COVID- 19) was first noted in December 2019 in Wuhan City (Hubei, China). The viral infection quickly spread worldwide, eventually causing a global pandemic. Following a detailed study of its biological properties, the virus was found to be of zoonotic origin and consists of a single-stranded ribonucleic acid (RNA) genome with a strong capsid. Based on this survey, it was concluded that the virus belongs to the coronaviridae family and was subsequently named 2019-novel coronavirus (nCOV). A person infected with 2019-nCoV may have no symptoms or develop mild symptoms, including sore throat, dry cough, and fever. If the human body hosts the 2019-nCoV for a long period, the virus can cause severe respiratory illness and, in the worst case, it can lead to death. There are four stages that are used to assess the virus's virulence in the human body. During the first four days of the infection, the patient is often asymptomatic. The second stage is the progressive stage which generally occurs between the fifth and eighth day following the infection, whereby the patient may develop mild symptoms. Stage three is known as the peak stage, which occurs between nine and thirteen days. The final stage is the absorption stage, whereby the load of the virus exponentially increases [1]. These observations were reported with clinical experimentation in Fig. 1 [2].

Figure 1:
An upsurge in the number of cases and death rate from January to July 2020 is depicted. The infection and death rate increased by approximately 10 5 within six months Due to the rapid surge in cases, healthcare systems are finding it increasingly difficult to cope with the demand and to provide timely vaccination [3]. This problem is being further exasperated by the shortage of medical supplies globally. In order to reduce the burden on healthcare systems, several preventive measures such as social distancing, proper sanitization, the mandatory wearing of masks in public places, and lockdowns have been implemented worldwide to reduce the spread. Despite the implementation of all these measures, the mortality rate from the disease is still high in various countries. According to the Chinese National Health Commission (NHC), as of February 4 th , 2020, the mortality rate from the disease was 2.1% in China and 0.2% outside of China. The mode of spread of the virus in asymptomatic cases remains controversial [4,5]. In order to identify COVID-19 in an asymptomatic person, precise and proper diagnostic tests are required. The diagnostic tests are typically performed by collecting samples from the individual patient for testing in a laboratory or at a point of care testing center [6]. Manual testing is time consumes and labor-intensive. Therefore this method is not suitable to obtain a fast diagnosis during a pandemic. Computed tomography (CT) and chest X-ray (CXR) can be used to detect and assess the severity of the lung damage caused by the viral infection. However, a radiologist needs to analyze these images manually, which is time-consuming. Artificial intelligence (AI) can be used to develop algorithms to automatically assess the lung damage caused by the virus [2,7]. The findings for the COVID-19 infection in CXR or chest CT vary from person to person. However, two common hallmark imaging features observed in infected patients were bilateral and peripheral ground-glass opacities and peripheral lesions with a rounded morphology [2]. These distinct features facilitate the use of machine vision learning models to automatically detect COVID-19 lesions on either CXR or CT images. However, traditional methods do not preserve the contextual information of CT scan images. In view of this, this study aimed to develop a robust diagnostic model for COVID-19 detection on CXR images. The objectives of this study were to: • analyze the behavior and performance of various vision models ranging from inception to Neural Architecture Search (NAS) networks followed by appropriate model finetuning, • visually assess the behavior of these models by plotting class activation maps (CAMs) for individual networks, • determine the classification performance of the model by calculating the area under the curve (AUC) of a receiver operator curve (ROC), • improve the generalization of the model by combining the finetuned model deep learning with the shocked model (stacked ensembles technique).

Previous Works
Numerous studies evaluated the use of deep learning methods for the automatic detection, classification, feature extraction, and segmentation for COVID-19 diagnosis from CXR and CT images. This study discusses the relevant applications of pre-trained deep neural networks that prompt the key aspects to impact COVID-19 detection and classification. Fan et al. [8] proposed the use of the deep learning network Inf-Net for the segmentation of COVID-19 lesions on transverse CT scan images. This network architecture utilized Res2Net as a backbone and obtained a dice score of 0.682. A similar semi-Inf-Net model attained a higher dice score of 0.739. Oh et al. [9] implemented two different approaches, global patch matching, and local patch matching, for segmentation and classification. Their method used ResNet-18 as the backbone to classify four different types of lung infections similar to that of COVID-19. Their algorithm obtained an accuracy score of 88.9% and specificity of 0.946 on randomly cropped patches using a local approach. Rahimzadeh et al. [10] constructed the 8-phase training concatenating Xception and ResNet-50 architectures. In each phase, samples were trained using a proper stratification to overcome class imbalance for 100 epochs. This model attained an overall accuracy score of 91.4% by five-fold cross-validation. Ozturk et al. [11] proposed a Dark-CovidNet model for binary and tri-class classification of CXR images infected with COVID-19. This model was trained by constructing a deep neural architecture with a series of convolutional layers and max-pooling layers. This method attained accuracy scores of 98.3% for the binary classification and 87.2% for the tri-class classification on five-fold cross-validation. Apostolopoulos et al. [12] applied transfer learning using diverse pre-trained architectures on two different datasets for the classification of COVID-19 CXR images. Their transfer learning methodology attained an accuracy score of 98.75% using VGG-19 pre-trained weights for binary classification and an accuracy of 94.7% for the MobileNet-V2 CXR images classification consisting of three classes. Li et al. [13] proposed the CovNet network by training a deep learning model with ResNet-50 as a backbone for sharing weights and attained an accuracy of 96%. Khan et al. [14] designed the CoroNet-architecture with Xception as an underlying weight-sharing model. This model achieved an accuracy score of 99% through binary classification, 95% when using three non-identical classes (one class belonging to , and 89.6% for four variant classes following a four-fold cross-validation framework. Wang et al. [15] proposed the use of COPEL-Net to segment COVID-19 pneumonia lesions from CT images. The novel dice loss combined with a MAEloss for generalization was used to reduce noise and minimize the foreground and background imbalance for the segmentation task. This diagnostic frame obtained a dice score of 80.72 ± 9.96. Most COVID-19 classification and segmentation on CXR and CT images described in the literature are based on deep neural networks. The advantage of deep neural networks is that they provide a versatile weight-sharing mechanism, thus improving the performance of the algorithm. Therefore, this study aimed to develop a robust diagnostic COVID-19 model using CXR images. The objectives of the study were to: • examine the behavior and efficiency of different deep learning vision models ranging from Inception to NAS networks, using the proper finetuning procedure, • visually assess the behavior of these models by plotting class activation maps (CAMs) for individual networks,

Dataset Description
A total of 2905 CXRs were obtained from various databases, including the Italian Society of Medical Radiology (SIRM), ScienceDirect, The New England Journal of Medicine (NEJM), Radiological Society of North America (RSNA), Radiopaedia, Springer, Wiley, Medrxiv, and other sources Fig. 2. The complete source list of the COVID-19 CXR image samples is available in the metadata file [16]. These images were reviewed by an expert radiologist. Eight percent (n = 219) of these images were from patients infected with COVID-19%, 46% (n = 1341) of the images were from healthy persons, and the rest of the images were from patients suffering from either bacterial or viral pneumonia (n = 1345) [16]. The data was then divided into 75% training (D train ) and 25% testing datasets (D test) . Due to the small number of CXRs with COVID-19 lesions, stratified random sampling was used to ensure that all three diagnoses were equally represented in both training and testing datasets and hence minimize the risk of introducing class imbalance in the data distribution.

Convolutional Neural Networks for Feature Extraction
Convolutional neural networks (CNNs) are increasingly being used in computer vision to detect, classify, localize, and segment normal and pathological features from medical images [17]. The use of CNN increased widely following its application in large-scale image recognition challenges (ILSVRC-2010). In this challenge, AlexNet [18] made use of a deep CNN and resulted in the lowest detection error rate. This motivated researchers to make use of this technology to develop multidisciplinary high-end applications [19]. The CNN architecture can be modified significantly by manipulating the width, depth, and channels (activation-maps) to further improve the performance of the model with appropriate generalizations. Furthermore, the model's performance can be further improved through the manipulation of parametric weight sharing from one network to another network. This technique facilitates the feature extraction procedure in most networks, eventually reducing the computational and training cost [20,21]. Following the successful implementation of AlexNet, numerous other CNNs were developed. In the following section, the advantages and limitations of each CNN are discussed.

Inception
The Inception architecture is designed with a novel ideology module. This network architecture is trained by widening layers to increase the depth of the network depth with a few computational parameters. There are two versions of the architecture, including a naive and a dimensionality-reduced. The Inception module consists of three levels. The bottom levels of inception feed into four different layers stacked by width. The intermediate layers extract spatial information individually and correlate with each layer. The top layer concatenates all the intermediate layer's feature maps to maintain a hierarchy of features to improve the perceived performance of the network [22].

VGG-Nets
After Inception, VGG networks were developed by a sequential convolutional layer with a pooling layer. The sequential depth of the models ranged from 11 to 19 layers. The appropriate use of the max-pooling layers in 16 and 19 layered VGG-Nets is essential for spatial sub-sampling and the extraction of generic features at the rearmost layers. VGG-Nets use small receptive fields of 5x5 and 3x3 to capture small features, eventually improving their detection precision accurately. The generalizability of the model for highly correlated inputs can be further improved by finetuning the learning application schedules to decrease the learning rate [23].

Res-Nets
The Res-Nets were developed to address the problem of vanishing gradients by imparting identity mapping in large-scale networks. They reformulated deep layers by aggregating learned activations from a prior layer to form a residual connection. This residual learning minimizes the problem of degrading and exploding gradients in the deeper networks. These residual connections help in addressing learned activations from preceding layers, maintaining a constant information flow throughout the network, and eventually reduce the computational cost [24][25][26].

Inception-Res-Nets
This network was inspired by the Inception network modules and identity mappings from ResNets. This method integrates dimensionality-reduced Inception modules with sequential residual connections hence increasing the learning capability of the network while reducing its computational cost. This provides better generalization ability when compared to various versions of the ResNet and Inception Networks [25].

Xception
This network was proposed to compete with the Inception network to reduce its flaws. The simultaneous mapping of spatial and cross-channel correlations guides allows for improved learning with small receptive fields and improves perceptive ability. The depth-wise separable convolutional layers enhance the learning through detailed feature extraction. These networks are computationally less expensive and perform better than the Inception network [27].

Dense-Nets
These densely connected CNNs are motivated by the residual connection of Res-Nets and imposed long-chained residual connections to form dense blocks. In Dense-Nets, for N layers, there are N(N+1)/2 connections (including residual connections) that enhance the network's capability for extracting detailed features while reducing image degradation. The sequential dense and transition blocks provide a collection of knowledge, and a bottleneck receptive field of 3x3, eventually improving its computational efficiency. The finetuning of larger weights improves generalization in deeper networks with a depth ranging from 121 to 201 layers [28].

Mobile-Nets
Mobile-Nets were designed for mobile applications under a constrained environment. The main advantage of this network is the combination of inverted residual layers with linear bottlenecks. The constructed deep-network accrues a low-dimensional input, which eventually expands by elevating dimensional space. These elevated features are filtered via depth-wise separable CNNs and are further projected back onto a low-dimensional space using linear CNNs. This contribution reduces the need to access the main mobile application memory, thus providing faster executions through the use of a cache-memory [29].

Nas-Nets
Nas-Nets make use of convolution cells by learning from distinct classification tasks. The design of this network is based on a reduced depth-wise stacking of normal cells, hence providing an appropriate search space by decoupling a sophisticated architectural design. This adaptability of Nas-Nets enables it to perform well even on mobile applications. The computational cost is significantly reduced, and its performance can be improved by enhancing the depth [30].

Deep Stacked Ensemble Method
This deep-stacked ensemble method was evaluated by classifying COVID-19 database inputs into a tri-class and a binary class, as shown in Fig. 3. Various samples were first considered and pre-processed to a specific resolution of 224×224×3 of the COVID-19 dataset. These pre-processed images were then fed into a variety of deep networks that use different paradigms to extract features from latent dimensions. The extracted feature vectors are then evaluated, and the two best-performing models are selected to form a stacked ensemble. The COVID-19 class is given more weight in this ensemble, which was assessed by classifying the feedback into a tri-class and a binary class.

Finetuning of Neural Networks
Deep learning algorithms can accurately detect pathology from bio-medical imaging to human-level precision. The CNNs provide numerous advantages for feature detection in medical imaging. There are two methods that can be used to design neural architectures for medical imaging. The first method involves designing a novel architecture by overhauling loops in existing architectures by training it end-to-end. The second method involves model finetuning by either transferring the weights of a pre-trained model (transfer of weights) or by retraining an existing pre-trained architecture.
The training of an end-to-end CNN requires proper initializations, which can be computationally expensive. On the other hand, the transfer of weights from the pre-trained models for a similar problem statement can be useful to reduce the computational cost. However, they may not extract the invariances if the class samples in the problem statement are not trained at least once. For example, a pre-trained network on Imagenet may not be able to extract the invariances in CXRs if these samples are never seen or trained. This means that the model may end up capturing unwanted features on the CXR, leading to an inaccurate classification. In order to overcome this problem, the model is fine-tuned to obtain the appropriate features. Fine-tuning of the model is extremely important in medical imaging when the sample size is small, leading to class imbalance [31]. Hence, the existing models, starting from VGG-Nets to Dense-Nets, were all finetuned to extract invariant features and discriminate the COVID-19 class from the remaining. The fine-tuning for individual models were performed as per Algorithm-1. The major parameters considered for finetuning in our methodology were learning schedules and batch sizes. The algorithms were finetuned by constricting the noise caused during the training process to reduce the risk of misleading the model if not trained with appropriate initializations.
The D train and D test samples were inserted into each model to capture latent feature vectors. A feedforward neural network was built to classify the extracted feature vectors, and all models were fine-tuned using Algorithm 1. The final extracted feature vector consisted of different threedimensional shapes according to the model. These latent representations were then classified by attaching a dense layer consisting of 256 neurons followed by dropout [32] and batch normalization [33] of the layers for regularization. The final layer consisted of a softmax activation layer with "c" neurons, whereby "c" represents the number of classes. The dropout percentage was set to 30%. A generalization assessment was performed for all individual models. ReLU [34] was used for the non-linearity construction of the model architecture for all the layers except the final layer, whereby feed was forwarded by softmax. Glorot-normal was used for the initializations of most of the layers [35]. The initializations with appropriate activations resulted in the extraction of the following intricate, deep feature layers.
All models were carefully finetuned, and their performance was evaluated using various performance metrics. The generalizations provided by the finetuned models are summarized in Tab. 1. All the models performed well and had a similar overall performance Tab. 2. The classwise performance of the model is also summarized in Tab. 2. The classes are coded as C-0, C-1, and C-2, indicating COVID-19, normal, and pneumonia, respectively.
In the design of medical diagnostic prediction models, receiver operating characteristics (ROC) analysis is essential for analyzing the model performance. The area under the curve (AUC) of the ROC of a classifier determines the diagnostic stability of the model. This AUC-ROC curve is insensitive to the alterations in the individual class distributions [36]. A ROC curve for each model was therefore plotted, as shown in Fig. 4. The feature extraction ability of models varied widely, as not all models were capable of recognizing features pertaining to COVID-19 lesions.
A prediction model for medical imaging needs to have a high sensitivity and specificity. A clinically useful COVID-19 model based on CXR needs to be able to differentiate between COVID-19 from other infections. However, the distinction of CXR lesions caused by COVID-19 as opposed to other infections can be quite challenging. CAMs were therefore applied to all CXR input images [37]. CAMs apply global average pooling for bottleneck activations in CNNs and provide a visual understanding of discriminative image regions and/or the region of interest. CAMs provide a visual illustration through the use of heat maps of the features extracted by the models to make predictions. Therefore CAMs provide a clear understanding of whether the acquired features are distinctive of a COVID-19 lesion, as illustrated in Fig. 5.
The CAMs analysis shows that some of the models extract the peripheral and bilateral ground-glass opacities while some of the other models also extracted the rounded morphology typical of COVID-19 lesions [38]. Since both features were deemed essential for an accurate diagnosis, the models that provided the highest generalization and extracted different features according to the CAMs analysis were used to develop the neural model averaging or neural stacked ensembles models.

Model Averaging
Model averaging is the process of averaging the outcomes of a group of networks trained on a similar task or the same model trained on different parameters. The model averaging improves the generalization of the models by aggregating their predictions. The generalization for the model was obtained by minimizing the loss during stochastic optimization using equation Eq. (1), whereby x and y are features and ground truth class labels of particular data distribution. If f n is an n th neural architecture that predicts the class label for a given feature set (where, n = 1, 2 . . . N), the mean squared error for the loss function can be minimized as follows: where f N represents the final neural architecture.
Similarly, weights can be assigned to individual models based on their prediction performance. These weights are then applied to the appropriate models to obtain aggregated generalization. This is known as weighted model averaging. In the case of model averaging, the models are equally treated by assigning the individual performance of the model to each network. This means that the weighted model averaging provides importance to the required models and discards the poorly performing models.
The generalization provided by the committee of the neural models improves when compared to that of model averaging and weighted model averaging. Hence, the models were stacked to improve the generalization ability of the model.

Stacked Ensembles
The stacked ensemble integrates or groups different models to provide aggregated generalization by mapping the output predictions onto a logit function. Instead of averaging the weights to the grouped models, logistic regression or multi-class logit was applied to map the predictions. Therefore, the predictions were gathered, and a logistic regression was applied to them or built at the end-to-end neural model that applies softmax non-linearity as final activation [39,40]. The generalization improvements provided by the stacked ensemble (using neural networks) were mathematically described as follows.
Our network was first considered to be a function that predicts a certain input x, where our true function is T i (x) and approximated function is f (x) ∀i = 1, 2, 3, . . .n. Suppose, where r i is the generalization error ∀ i = 1, 2, 3 . . .n , whereby n represents the number of neural networks to an ensemble. where So, the average individual error settled from the networks can be estimated as follows: The ensemble learning of the grouping variant networks is presented in the following equations: Estimated error resided by stacking these ensembles: Suppose, From Eqs. (13) and (8) If individual networks did not correlate themselves, the stacking ensemble was reduced by n factors using the original generalization attained from the individual networks. However, for most scenarios, a correlation in generalization occurs, leading to an increase in the generalization error to a certain extent.
To understand this scenario, r ij = 0 is considered. So, This constant 'ε' is an additional error caused due to covariances underlying the perception of individual networks.
With this knowledge, it is clear that the stacked ensembles can outperform single networks in terms of generalization. As a result, six different neural network committees were formed by multiplying the number of neural networks described in Tab. 3, which ranged from 2 to 13 networks (all). These ensemble networks were evaluated using the standard classification metrics. A small neural architecture was attached to the committee of networks to adjoin the connected layer fully. This fully CNN consists of 16 neurons with a dropout of 30% for regularization. The final activations were pushed with softmax non-linearity, which consisted of three neurons describing the class predictions pertaining to each individual class. The results obtained by the proposed stacked ensembles are described in the next section in detail. Stacking all the models in Table. 1

Results
As mentioned, the generalization error obtained by a committee of neural networks is always less to that of a single neural network. Six variant committees of networks were selected and combined as described in Tab. 3. The classwise classification metrics utilized to understand the behavior of a specific COVID-19 class are illustrated in Tab. 4.
A comparative study was then performed to compare the performance of the proposed network with other existing models described in the literature Tab. 5.
Our designed generic training algorithm facilitates the training process by acquiring faster convergence and with low computations (Iterations). During the training process, the batch size and learning rate are increased cautiously for each iteration to obtain a balanced criterion, as explained by Smith et al. [41]. As noise during training can be reduced by properly choosing the batching parameters fed into the network, the learning rate and momentum of the optimizer were assigned a faster search.   The noise due to training is theoretically represented as follows in Eq. (18).
Here we assumed a constant momentum. A training algorithm was developed to conceptualize the noise constraint. Although a decaying learning rate can decrease the noise, it gradually increases the computational time for training. On the other hand, lowering the batch size can also reduce noise but comes at the cost of lowering the generalizing capacity of the model. These problems were overcome by developing an algorithm that increased the batch size during the specified iteration and cautiously increasing the learning rate as follows. The algorithm was first iterated for 16278 steps (iteration 1), whereby the learning rate was set to 10 −4 by sending 15 samples as a batch at a time. In the next iteration (iteration 2), the batch size was increased by 50%, and the learning rate was increased tenfold. In order to maintain a consistent trade-off between generalization and faster convergence, the batch size was increased by a factor of 150% (to that of initial), and the learning rate was tuned as per the preceding iteration. During the experimentation, it was found that the proposed training procedure led to a faster convergence by training using only a few steps (approximately 20 epochs).
Appropriate training with fine-tuning of the ensemble was therefore critical to obtain these insightful outcomes. The final Ensemble-6 model had the highest performance when compared with the other method, with an accuracy score of 99.175%. Ensemble-1 and Ensemble-2 attained an accuracy of 98.487% and 98.762%, respectively. When taking into consideration only the COVID-19 class, the precision rate was at least 97.674%, but the recall rate was lower. The highest and lowest recall rates were 89.795% and 69.387% and were obtained by Ensemble-6 and Ensemble-4, respectively. However, due to the small sample of the COVID-19 class in our study, it was difficult to extract additional invariant features to improve the performance of the model further.

Limitations
In this study, we observed that the stacked ensemble was slightly inefficient when a poorperforming model was included. The DenseNet-201 model evaluations were not always finetuned correctly, and the network depth was not always appropriate, leading to a high generalization error. The COVID-19 results were not always included in the single model based on the features derived from the individual models. The ensemble method offers more generalization, but the combination of multiple models increased the computational cost, which is unnecessary for smallscale computational systems (such as Ensemble-6). As a result, in real-world scenarios, small, quick, and efficient models such as Ensemble-1 and Ensemble-2 are advantageous. The progression of the virus can be visualized better on Chest CT axial images. However, there is a chance of missing disease progression on CXR [38], which could be dangerous. Therefore future studies should focus on the development of models that can predict disease progression on CXR.

Conclusion
In this study, various COVID-19 classification models were evaluated and compared using different classification metrics. Furthermore, a learning framework for finetuning these models was proposed, and their bottleneck activations were visualized using CAMs. The AUC-ROC curves were closely examined, and the output of each class was illustrated visually. These finetuned models were then stacked to outperform previous models and include a broad range of generalizations. The ensemble models achieve an accuracy score of 97.66 percent in the worstcase scenario. Even after finetuning for class imbalance, the models were found to have a high generalization ability. The least error rate obtained by the outperforming model, built by stacking all the finetuned models, was 0.83%. The stacked ensembles method improved the performance of the model and could therefore be used to improve the prediction accuracy of the diagnostic models in medical imaging.