Robust System Identification: Finite-sample Guarantees and Connection to Regularization

${\color{black}\text{1. Theoretical results}}$

• We agree with the reviewer that our analysis was limited to the linear case. As we acknowledged in the paper, despite its apparent simplicity, the linear system remains challenging to analyze. Most research in statistical learning for system identification has focused on linear systems, as this setting allows for more tractable theoretical analysis.

• However, we’d like to remark that a nonlinear setting, in fact, spans a very broad regime. Several studies have provided non-asymptotic guarantees for subclasses of nonlinear systems that exhibit near-linear behavior. While we appreciate and could have adopted these approaches, their analyses fundamentally rely on techniques similar to those used in linear settings, leveraging near-linear behaviors. Our priority in this paper, however, was to convey the main insight that “regularization or robustification improves empirical performance with nearly no additional theoretical cost (i.e., nearly the same theoretical error)”, encouraging practitioners to adopt our approach over the standard LSE.

${\color{black}\text{2. Comparison with other Schatten norms AND the Wasserstein robust models}}$

We thank the reviewer for the excellent question.

• As mentioned briefly in the proof of Proposition 1, the choice of $p$ in the Schatten $p$-norm—or equivalently, $q$ in our SDP formulation (7)—has only a minimal impact on the theoretical guarantee, due to the "equivalence of norms."

• However, in practice, the Schatten norm parameter $q$ in our SDP formulation (7) can be selected based on how the eigenvalues of the underlying system are distributed. To illustrate this relationship more clearly, we conducted stylized experiments: ${\color{maroon}\text{for the detail of the experiment, please check Fig. 1 in } **Supplementary Material**} \text{ or Appendix A.8 in the updated manuscript.}$

• We conducted adaptive control tasks with different choices of $q$: ${\color{maroon}\text{for the detail of the experiment, please check Fig. 2 in } **Supplementary Material**} \text{ or Appendix A.8 in the updated manuscript.}$

• We further conducted the following experiment where we compared ours with the SOTA Wasserstein model: ${\color{maroon}\text{for the detail of the experiment please check Fig. 3 in } **Supplementary Material**} \text{ or Appendix A.6 in the updated manuscript}$.

${\color{black}\text{3. The optimal solution to the true problem}}$

• We believe that the validity of Eq. (5) is a straightforward extension of the consistency of OLS, which is well-documented in many textbooks (e.g., [Shao, 2008]).
• The following confirms that Eq. (5) holds for $T=1$. For the general case where $T>1$, the same idea applies with only more cumbersome recursion. So we omit that in the rebuttal; a similar derivation for a linear case can be found in [Dean et al., 2020].

First of all, the minimization problem in Eq. (5) is equivalent to

which we call the true LSE problem. Recall that our linear system is

For the case of $T=1$, the objective function is $f(\theta)=\mathbb{E}[\Vert {x}\_1+\theta \phi({z}\_0)\Vert_2^2]$. Then, we have

$\hspace{3cm} \theta^{\prime} \in \underset{\theta}{\operatorname{argmin}} f(\theta)$ and $\frac{\partial}{\partial \theta} f(\theta)=0 \Rightarrow 2 \theta^{\prime}\mathbb{E}[ \phi\left({z}\_0\right) \phi\left({z}\_0\right)^{\top} ] = 2 \mathbb{E}[ {x}\_1 \phi\left({z}\_0\right)^{\top} ].$

$\theta^{\prime}= \mathbb{E}[ {x}\_1 \phi\left({z}\_0\right)^{\top} ] \left( \mathbb{E}[ \phi\left({z}\_0\right) \phi\left({z}\_0\right)^{\top}] \right)\^{-1}$. Plugging $x_1=\theta^{\star} \phi\left(z_0\right)+w_0$, we have

$\begin{aligned} \theta^{\prime} & = \theta^{\star} + \mathbb{E}\left[ w_0 \phi\left({z}\_0\right)^{\top}\right] =\theta^{\star}\end{aligned}$

{\color{black}\text{4. Definition of \Sigma_w}}

• Thanks for your careful review! $\Sigma_w$ was meant to be the covariance of the noise term. We will revise the paper accordingly.

${\color{blue}\text{References}}$

• Shao, Jun. Mathematical statistics. Springer Science & Business Media, 2008.

• Dean, Sarah, et al. "On the sample complexity of the linear quadratic regulator." Foundations of Computational Mathematics 20.4 (2020): 633-679.

${\color{black}\text{1. Motivating examples of the single trajectory setup}}$

We thank the reviewer for the excellent question.

• We believe that the single-trajectory setting more accurately reflects the realistic data collection process: in dynamic systems, data collection naturally results in non-iid data since each new observation depends on the current state.

• In this context, the iid data assumption is more restrictive, namely, the iid assumption is that we have $N$ iid trajectories, $\\{x_t^i\\}\_{t=0}^T$ for $i = 1, \ldots, N$. As noted in [Dean et al., 2020], the iid assumption implies that we somehow have the ability to reset the “unknown” system to a state independent of past observations (e.g., $x_0 = 0$).

• While resetting a system may be straightforward in certain robotics experiments—where an agent can easily be returned to an initial state such as a specific position or angle—this is not always feasible. For example, in safety critical missions such as drones deployed in space or underwater, data-collection is very limited.

• This is also related to the advantage of conventional pipelines. A model-free approach can provide very flexible policies that a model-based approach may not come up with. However, the model-based approach is more data-efficient, hence, more suitable under a low data regime.

${\color{black}\text{2. Comparison with the SOTA method AND in a large scale system}}$

We really appreciate the constructive feedback from the reviewer.

• We’d like to point out that wind prediction is a real-world problem that has been considered by various works.

• Following the suggestion we conducted the following experiments with a large scale system with dimensions being $50\times100$: ${\color{maroon}\text{please check Fig. 5 in } **Supplementary Material**} \text{ or Appendix A.7 in the updated manuscript}$. The results not only align but also corroborate our previous experiments.

• We further conducted the following experiment where we implemented the SOTA Wasserstein model: ${\color{maroon}\text{please check Fig. 3 in } **Supplementary Material**} \text{ or Appendix A.7 in the updated manuscript}$.

$\quad$ - Our method achieves better performance in mean norm error. We also demonstrate the computational advantage of our model: the run time for our SDP formulation remains constant with sample size, whereas the run time for the Wasserstein model scales with sample size.

${\color{black}\text{3. Assumptions in Sec. 4}}$

• First, we’d like to point out that A2 is a standard assumption in this stream of research as described in the review paper [Tsiamis et al., 2023].

• That said, we agree that A6 may be too strong. In fact, A6 is a simplifying assumption intended to improve the readability of our paper and avoid overly cumbersome notation in our theoretical derivations. A more realistic assumption in this line of research is the following: we have access to some policy $K_0$ that stabilizes the unknown system. As stated in [Dean et al., 2020], in many cases a stabilizing policy can be found efficiently. With this assumption replacing A6, we can inject exploration noise into the stabilizing policy, i.e., $u_t = K_0 x_t + \eta_t$, where $\eta_t$ is sub-Gaussian noise, and derive similar guarantees to Proposition 1.

${\color{black}\text{4. Extension to unobservable state setting}}$

We thank the reviewer for the interesting question.

• As briefly mentioned in Conclusion, extending the proposed approach to unobservable systems is a direction for future research.

• In [Oymak and Necmiye, 2019], the authors consider an unobservable linear time-invariant system:

They show that the standard LSE can be used to recover the Markov parameters of the linear system and reconstruct $A, B, C, D$ from the Hankel matrix via the Ho-Kalman algorithm. Using this idea, we can derive a similar SDP formulation as in Eq (7) for such a system.

• While the extension to the unobservable setup is direct from an algorithmic development standpoint, we are unsure whether a new approach is required to derive theoretical results. It would be interesting to explore if we can establish theoretical guarantees for our approach in this setup.

${\color{blue}\text{References}}$

• Dean, Sarah, et al. "On the sample complexity of the linear quadratic regulator." Foundations of Computational Mathematics 20.4 (2020): 633-679.

• Tsiamis, Anastasios, et al. "Statistical learning theory for control: A finite-sample perspective." IEEE Control Systems Magazine 43.6 (2023): 67-97.

• Oymak, Samet, and Necmiye Ozay. "Non-asymptotic identification of LTI systems from a single trajectory." 2019 American control conference (ACC). IEEE, 2019.

We thank X6xt for the valuable feedback and Wdgh for providing further explanations and references. The rebuttal process has been very helpful for us to identify areas we may have overlooked. In particular, we are less familiar with the literature on improving the computational aspects of DRO using GPUs, so it’s great to learn about such efforts. Our comparison was performed assuming that the system operator has access only to standard off-the-shelf solvers.

• First, we’d like to clarify that the Wasserstein model used in our experiments is not an LP, but an SDP, since the objective function is squared errors. As the reviewer mentioned, we could implement the LP version of the Wasserstein model by changing the objective function to absolute errors. However, our experiments were intended to compare the Wasserstein "LSE" model with our robust LSE model. In this case, the Wasserstein model shares the same complexity as our formulation, even when the sample size is fixed. Unfortunately, we are not aware of Wasserstein LSE models that have LP or SOCP formulations.

• So, we certainly agree with the reviewer’s concern. In our tests with linear systems with a state dimension being 50, solving the SDP formulation using MOSEK (solver) took an average of 40 seconds. Therefore, for such large systems, we agree that an SDP formulation is not a practical choice to begin with.

• That said, as mentioned in our paper, the SDP formulation in Equation (7) can be reduced to an unconstrained quadratic program (QP) by setting $q = 1$. While this approach loses the flexibility of selecting the Schatten norm parameter, it certainly offers a reasonable tradeoff given the computational concerns. On the other hand, the Wasserstein model does not admit such regularization interpretation for a least squares loss function so its formulation remains not scalable. We will discuss the limitations of our approach and SDP formulations in the paper.

{\color{black}\text{1. The practical or theoretical consequences of varying Schatten norm parameter p }}

We sincerely thank the reviewer for the constructive feedback.

• As mentioned briefly in the proof of Proposition 1, the choice of $p$ in the Schatten $p$-norm—or equivalently, $q$ in our SDP formulation (7)—has only a minimal impact on the theoretical guarantee, due to the "equivalence of norms."

• However, in practice, the Schatten norm parameter $q$ in our SDP formulation (7) can be selected based on how the eigenvalues of the underlying system are distributed. To illustrate this relationship more clearly, we conducted stylized experiments: ${\color{maroon}\text{please check Fig. 1 in } **Supplementary Material**} \text{ or Appendix A.7 in the updated manuscript.}$

$\quad$ - We generate two $3 \times 3$ linear systems for the experiment. The system in the left subplot has uneven eigenvalues ($\lambda_1 = 0.9$ and $\lambda_2 = \lambda_3 = 0.1$), while the system in the right subplot has even eigenvalues ($\lambda_1 = \lambda_2 = \lambda_3 = 0.9$). We conduct 300 simulations across different choices of the Schatten norm parameter $q\in\\{1,2,\infty\\}$ and the radius parameter $\epsilon$ (on the x-axis), plotting the mean norm error (on the y-axis) to identify the optimal radius for each value of $q$. As shown, $q = 1$ performs the best for the system with uneven eigenvalues, while $q = \infty$ performs the best for the system with even eigenvalues.

• We further conducted adaptive control tasks with different choices of $q$: ${\color{maroon}\text{please check Fig. 2 in } **Supplementary Material**} \text{ or Appendix A.7 in the updated manuscript.}$

$\quad$ - For example, $**i)**$ Boeing 747 represents a system with uneven eigenvalues, whereas $**iii)**$ UAV corresponds to a system with even eigenvalues. We conducted 500 simulations to evaluate adaptive control tasks using different values of $q$ in Robust STABL (R-STABL). The results align with our interpretation in Fig. 1.

${\color{black}\text{3. Avoiding the curse of dimensionality}}$

We thank the reviewer for the clarifying question.

• We may mislead the reviewer in using the term “curse-of-dimensionality” when comparing our result with the Wasserstein robust model. As mentioned, the Wasserstein model encounters the curse-of-dimensionality because its error rate is $\mathcal{O}(1/T^{\frac{2}{n}})$ where $T$ is the number of samples and $n$ is the dimension of state space. We used the term “the curse-of-dimensionality” to describe this exponential dependence of $T$ (i.e., # of sample) on $n$, which is not present in the guarantee provided by our robust model. Addressing this exponential dependence is an open research topic that the robust optimization community tries to address. Although a different approach was proposed by [Gao, 2023] under the assumption of iid data and Lipschitz continuity with respect to $\theta$, this approach is not applicable to our system identification problem.

${\color{black}\text{4. The difference between the curse of dimensionality and the unavailable dependence of the system dimension in the burn-in}}$

We appreciate the reviewer for the careful feedback.

• As another reviewer suggested, $T$ scales with other problem-dependent parameters, most notably the dimension of the system, i.e., $n$ and $m$. We acknowledge that it would have been more informative to present the explicit upper bound for the radius parameter $\epsilon(\delta)$ in Eq. (11). We have simplified the explicit bound as shown below:

Here, the constants $c_1,\\;c_2,\\;c_3,\\;c_4,\\;c_5$ are universal constants.

${\color{blue}\text{References}}$

• Gao, Rui. "Finite-sample guarantees for Wasserstein distributionally robust optimization: Breaking the curse of dimensionality." Operations Research 71.6 (2023): 2291-2306.

${\color{black}\text{1. Comparison with the SOTA method}}$

We really appreciate the constructive feedback from the reviewer.

• Indeed, other reviewers suggested comparing our method against the fine-tuned SOTA Wasserstein robust model: ${\color{maroon}\text{please check Fig. 3 in } **Supplementary Material**} \text{ or Appendix A.7 in the updated manuscript}$.

$\quad$ - Our method achieves better performance in mean norm error. We also demonstrate the computational advantage of our model: the run time for our SDP formulation remains constant with sample size, whereas the run time for the Wasserstein model scales with sample size.

${\color{black}\text{2. Comparison with a marginally stable system}}$

Thanks for the excellent suggestion.

• Following the reviewer's suggestion, we further conducted the experiment with a marginally stable system and simulated the stability: ${\color{maroon}\text{please check Fig. 4 in } **Supplementary Material**} \text{ or Appendix A.7 in the updated manuscript.}$

$\quad$ - As the reviewer suggested, the results seem to illustrate that the stabilization rate and system errors are related.

${\color{black}\text{3. Comparison to the DRO literature}}$

• The analysis in [Gao, 2023] requires Lipschitz continuity of the loss function and boundedness of the feasible set. Neither is satisfied in our setting, and it's unclear how the proofs can be extended to our setting.

{\color{black}\text{4. The valid upper bound on \epsilon(\delta)}}

We thank the reviewer for the careful feedback

• We acknowledge that it would have been more informative to present the explicit upper bound for the radius parameter $\epsilon(\delta)$ in Eq. (11). We have simplified the explicit bound as shown below:

Here, the constants $c_1,\\;c_2,\\;c_3,\\;c_4,\\;c_5$ are universal constants.

${\color{blue}\text{References}}$

• Gao, Rui. "Finite-sample guarantees for Wasserstein distributionally robust optimization: Breaking the curse of dimensionality." Operations Research 71.6 (2023): 2291-2306.

Thank you for the reply. The authors have addressed my main comments satisfactorily, and I have increased my score. I am particularly satisfied by the % stabilized experiment (Figure 10), which concretely shows the "safety" improvement from using the proposed robust estimator for sysID + certainty equivalent control over naive LSE, even in the face of (nominal) marginal stability.

I do agree with Reviewer q1qq that misspecification is a very important desideratum that is largely left open. However, I'm willing to view this as too large a can of worms for this paper, as it is not obvious to me the best way to quantify "misspecification" in a non-trivial, control-oriented way beyond a naive norm perturbation. Furthermore, the aforementioned stabilization experiment partially assuages this point for me, as one can view certainty equivalent control on the identified dynamics as a proxy for performance under misspecification.

${\color{black}\text{1. Model misspecification}}$

We really appreciate the enlightening feedback from the reviewer.

• Poor out-of-sample (i.e., test) performance caused by small data is one of the primary motivations behind robust optimization. Numerous studies in robust optimization demonstrate that robust solutions enhance out-of-sample performance [Bertismas et al., 2018]. That said, we acknowledge that the logical flow might have been clearer for a general audience if we had begun with the regularization formulation. We will update our manuscript accordingly.

• However, we would like to emphasize that the ideas in this paper originated from an attempt to apply robust optimization to system identification, without initially anticipating its equivalence to the Schatten $p$-norm regularization. As the reviewer noted, we did not provide theoretical implications of robustness against model misspecification while the wind prediction experiment was intended to empirically demonstrate this robustness: the neural network model is an approximation of the dynamic of wind speed. In response to the feedback, we have considered the following.

Recall that our linear system (1) is

and the uncertainty set (6) is

Instead of (6), consider the following set:

As shown above, instead of a set of the matrix $\Omega$, we can express our uncertainty set as a set of nonlinear functions $\psi(\cdot)$. This alternative expression implies that the uncertainty set constrains the functional space of $\psi(\cdot)$. Therefore, we can interpret our Robust LSE problem in the following:

In this context, our problem provides robustness against model misspecification in the nonlinear function $\phi(\cdot)$ if $\phi(\cdot)$ is indeed unknown like the wind speed prediction example in our paper.

${\color{black}\text{2. Novelty}}$

• As proposed, our primary goal in this paper is to address the poor empirical performance of LSE under data scarcity and/or model misspecification. We appreciate the reviewer’s suggestion that regularization could be a natural choice to address these issues. However, to our knowledge, no prior papers are providing a theoretical non-asymptotic guarantee for regularized LSE in the standard setup. In this regard, we'd like to emphasize that our results are new while we acknowledge that the derivations rely on standard convex optimization theory and statistical learning theory.

• As briefly noted in our paper, most regularization methods rely on heuristic choices for the regularization parameter, often without a solid theoretical foundation. For instance, the parameter is frequently set to a fixed small value or made independent of the number of samples. In contrast, our results suggest that the regularization parameter should be data-dependent to achieve optimal performance.

• Specifically, Theorem 1 and Corollary 1—which establish the equivalence between regularization and robustness—represent a significant contribution, as this equivalence allowed us to adopt a distinct (i.e., robust) framework to derive non-asymptotic guarantees for regularization, thereby advancing theoretical understanding in the system identification community.

• Our results also have implications for the robust optimization community, as our robust model addresses the curse of dimensionality, which remains a challenge for the state-of-the-art Wasserstein model (${\color{maroon}\text{please check Fig. 3 in } **Supplementary Material**} \text{ or Appendix A.6 in the updated manuscript}$ for the numerical experiment.) In this way, we believe that our theoretical contributions are relevant to both communities by bridging system identification and robust optimization.

${\color{blue}\text{References}}$

• Bertsimas, Dimitris, Shimrit Shtern, and Bradley Sturt. "Two-stage sample robust optimization." Operations Research 70.1 (2022): 624-640

Dear Reviewer q1qq

Thank you for your thoughtful feedback.

We would appreciate it if you could let us know whether our response addresses your concerns. If so, we kindly ask you to consider raising the rating of our work.

Thank you again for your time and effort! We are looking forward to discussing any additional concerns you may have.

Best wishes,

Authors