We thank the reviewer for the careful and thorough review and provide the explanation below:

• Trajectory length: The reviewer is absolutely correct that by picking a minimally satisfactory $T$ with $m \approx k$ and $T \approx 3m$, simulations have shown we can increase the probability of getting accurate system identification and stabilization policy almost surely by increasing $M$. This roughly matches the order in the theorem ($m$ scales linearly with $k$), though we do acknowledge that the theorem has very large numerical constants and $m$ needs to depend on other constants as well. We have added discussions on this in the paper.

• Bounded Hankel around line 67: The reviewer is right that we should have been more careful with word choice. While it is true that in a strictly finite observation scenario, Hankel matrix is always finite and bounded, estimating Hankel matrices accrately becomes much more difficult when the system is unstable. As the $(i,j)$-th entry of the Hankel matrix is $CA^{m+i+j-2} B$ for some $i,j \in \mathbb{N}$, if $A$ is stable, $\Vert CA^{m+i+j-2} B\Vert$ decreases with increase in $i,j$, so the error of each entry decreases as the trajectory gets longer. However, if $A$ is unstable, $\Vert CA^{m+i+j-2} B\Vert$ increases, so any error in system estimation of $A$ (and by extension $B,C$) would be amplified, making the problem more difficult. We will clarify this in the revision.

We hope the above have answered your questions and provided some more context. Please let us know if you have further questions or suggestions.

We thank the reviewer for the careful and thorough review and provide the explanation below:

• Scale of constants: The constant $\epsilon_*$ is a stability margin linking the necessary condition for the existence of a stabilizing controller, i.e. $\Vert K(z)F_K(z)\Vert_{H\_{\infty}} < \frac{1}{\Vert \Delta(z)\Vert_{H\_{\infty}}}$ as in Lemma 3.1 and the actual bound on $\Vert K(z)F_K(z)\Vert_{H\_{\infty}}$ in Assumption 5.2 which requires some margin on top of Lemma 3.1. If we interpret $\epsilon_*$ to be in the same order as $\Vert\Delta(z)\Vert_{H_{\infty}}$, then it does not inherently scale with $n$ as $\Vert \Delta(z)\Vert_{H_{\infty}}$ characterizes the input-output frequency response, related to the input-output dimension, not directly related to the internal state dimension. On the other hand, $L$ covers all other constants needed in Theorem 5.3, it upper bounds the matrix norm of all of the following matrices $\Vert A\Vert, \Vert B\Vert,\Vert C\Vert, \Vert N_1\Vert, \Vert R_1 B\Vert, \Vert CQ_1\Vert$. If we take $A$ to be from a large-scale networked LTI systems, i.e. with very large $n$ but $A$ admits a sparsity pattern to a sparse graph, then $\Vert A\Vert$ scales with the maximum degree of the graph (assuming each entry of $A$ is $O(1)$).

  The dependence on $n$ is also shown in simulation, i.e. as we increased the state dimension $n$ with randomly generated matrices, the length of trajectories $T$ and number of trajectories $M$ stayed the same.

  Finally, we point out that even if $\epsilon_*, L$ depends on $n$ polynomially, its effect on the final complexity bound is only $\log n$ as $\epsilon_*,L$ only shows up in $\log$ terms.

• $M,T$ tradeoff: Thanks for your suggestions. The $M$ and $T$ trade off is indeed the idea we want to analyze in the simulation. Although the preliminary result shows that fewer trajectories (smaller $M$) requires longer trajectories (larger $T$) and shorter trajectories (smaller $T$) requires more trajectories (larger $M$), we did only include one simulation for a single fixed $T$ and varying $M$ and another simulation for a single fixed $M$ and varying $T$. We did not include more results because it seems the trend of the plots look similar to what is present. However, the reviewer is right in pointing out that there might be value in collecting the percentage of successful stabilization for a fixed $n$ across different $M,T$ pairs. We are currently completing this simulation and will include that in the final version.

We hope the above have answered your questions and provided some more context. Please let us know if you have further questions or suggestions.

We thank the reviewer for the careful and thorough review and provide the explanation below:

• Regarding $\delta$: We thank the reviewer for pointing out the typo in the statement of Theorem 5.3. It should have been "with probability at least $1-\delta$", instead of "with probability at least $\delta$".

• Assumption 5.2 and no marginally stable eigenvalue: Both assumptions are necessary for this unstable+stable decomposition to work. As discussed at the end of Section 5, Assumption 5.2 comes from the Small Gain Theorem, which is an if and only if condition. If Assumption 5.2 is not true, it can be shown that no matter what control law $K(z)$ the user picks, there always exists some $\Delta(z)$ that would make the system unstable [53]. Regarding marginal stable eigenvalues, as our results depend on the separation of the unstable and stable components, we do not allow marginal eigenvalues. That being said, we do allow the eigenvalue to be arbitrarily close to $1$, and its impact on sample complexity is characterized in our sample complexity bounds, discussed under the ''Dependence on system theoretic parameters'' paragraph after Thm 5.3.

  Regarding your point ''The work doesn't give a full solution to the stabilization problem:'', we wish to point out that the learn-to-stabilize problem has been shown as a challenging problem with various exponential lower bounds proven in [9] [48] and further in paper ''Learning to Control Linear Systems can be Hard''. Our result essentially carved out a significant subset of problems (those satisfying assumptions in the paper) for which learning-to-stabilze can be easy. In fact, our assumptions is minimal in order to break the hardness of stabilization proposed in [9,48] if one only stabilizes the unstable subspace of the system. While we do acknowledge there is more to do towards ``full solution to the stabilization problem'', our results already make a significant step towards this problem as all prior results actually show learn-to-stabilize is hard.

• Process noise: The paper's theoretical guarantee works when there is no process noise, but observation noise is allowed (Equation (1), (14), (26) all include the observation noise $v_t \sim \mathcal{N}(0,\sigma_v)$). Process noise is only included in Equation (26) as $w_t$ to show that the proposed method works under process noise in simulation. As the reviewer points out, introducing process noise would enlarge the error of estimating $\mathcal{H}$ with an additional error term in Equation (31).

  We clarify that the main technical novelty of the paper is the stable + unstable decomposition (including techniques to separate the unstable component from the lifted Hankel), whereas the system ID part to get the Hankel is not the main focus. We chose to assume process noise is $0$ to keep the system ID part cleaner, avoiding complicating the paper with unnecessary details. If process noise are to be included, then as long as the Hankel estimation error can be bounded, the rest of the proof would work in a similar way. Here is the outline of what the proof would look like:

  If there is process noise $w_t$, a matrix $\Psi = [0, A, A^2, \dots,A^T]$ will have to appear in the proof, so the estimation error of $\hat{\tilde{H}}$ will contain a term that increases exponentially with $T$. This is to be expected, as process noise would be amplified by the unstable part of $A$. Once we bound the Hankel estimation, the rest of the proof will be exactly the same as in this paper. In order to finish proof in this direction, a fine balance of $T$ needs to be selected, so that it needs to minimize the current terms in (109), but also not make a term containing $\Psi$ blow-up, so more weight needs to be applied to the number of trajectories $M$, instead of $T$. However, this fine balance seems only to be present on the above theoretical analysis, as in simulation, we do include process noise we see that (a) the number of samples needed for stabilization does not increase with $n$, confirming our approach works under process noise; (b) a longer trajectory still helps with stabilization even with process noise, meaning the fine-balance regarding the choice of $T$ is not present in our simulation. We leave how to improve the theoretical analysis to incorporate process noise as a future direction of this paper.

• Initial condition: The paper does assume the system is reset to $0$ at the start of each trajectory, but this is merely for the simplicity of proof. For simplicity, we will set $D = 0$ again in this part of the discussion. If $x_0 \neq 0$, then $y_{t+1} = CA^t x_0 + \sum_{i=0}^{t}CA^{t-i}B u_i$, with the additional term $CA^t x_0$ known as the free-response term. We can pick some $p$ such that $k < p <T$ and define $U_p = [u_0,\dots, u_{T-p}], U_f = [u_p, \dots, U_T]$ and $Y_p = [y_0, \dots, y_{T-p}], Y_f = [y_p \dots, y_T]$. We can then project the future outputs orthogonally to the row space of $Y_p$ by $\Pi^{\perp} = I- Y_p^\top (Y_p Y_p^{\top})Y_p$. Then, we can use the observation $\tilde{Y}_f = Y_f \Pi^\perp$ for the proposed algorithm. This is a quite popular method used in conventional control methods, such as N4SID. For details on how and why it works, please refer to Section 3 of ''N4SID: Subspace Algorithms for the Identification of Combined Deterministic-Stochastic System'' by Peter Van Overschee and Bart De Moor. This is indeed an important issue; we have added more discussion on it in the paper.

• Inconsistencies in $m,p,q$ btw theory and simulation: Theorem 5.3 offers the order of $m,p,q$, not the exact value we should pick (the $\asymp$ means only ignoring numeric constant as defined at the start of Section 2 ). These expressions provide theoretical insights, showing how $m,p,q$ should depend on system size $n$, unstable subspace $k$, observability, and contrability of the system. Making an exact computation of the omitted numerical constants might offer conservative values, as throughout the proof, many terms are relaxed and merged with the leading order term, making the final order-wise dependence clean at the cost of a conservative numerical constant. For this reason, in experiments, we do not follow the exact theoretical values, but follow trial and error instead. This is common practice in many theoretical literature in control and optimization, where the parameters chosen in the Theorems are conservative, and simulations use much more relaxed parameters.

• Usefulness of bounds for $m$: We disagree that ``the choice of $m$ ... does not appear to be very useful''. The choice of $m$ eventually translates to the sample complexity requirements $(3m+2)M$, which answers fundamental questions of how the sample complexity of stabilization depends on system theoretic parameters like $n,k$, eigenvalue gap, controllability gramian etc. These bounds provide insights like some systems are inherently easier to stabilize than others. For example, we see it is harder to separate the stable system from the unstable system if there are marginally stable eigenvalues. This is why $m$ is proportional to $\frac{1}{1-|\lambda_{k+1}|}$ and $\frac{1}{|\lambda_{k}|-1}$. Moreover, the less observable or controllable a system is, the smaller $\sigma_{\min}(\mathcal{G}\_{\text{ob}})$ or $\sigma_{\min}(\mathcal{G}\_{\text{con}})$ is, respectively, and the larger $\kappa(\mathcal{G}\_{\text{ob}})$ or $\kappa(\mathcal{G}\_{\text{con}})$ is, respectively. Therefore, the less observable or controllable a system is, the harder it is to stabilize the system.

  We acknowledge that in practice, the theoretical values for $m$ cannot be exactly computed, as it depends on the parameter of the unknown system. However, we point out that as long as we plug in conservative estimates of those parameters into $m$ (e.g. lower bounds for $\sigma_{\min}(\mathcal{G}\_{\mathrm{ob}}),\sigma_{\min}(\mathcal{G}\_{\mathrm{con}})$), the theorem still works. Further, as in many theoretical papers, the selection of those parameters in experiments may differ from those in theory. In fact, in our experiments, the choice of $m$ is very forgiving. Unless $m$ is too small (so we did not converge to unstable subspace enough) or too close to $T$ (so we did not leave enough samples for system estimation), we generally get quite good system estimation. Throughout the simulation, we did not need to do much hand-picking on $m,p,q$. We have added more discussions on this issue in the paper.

We hope the above have answered your questions and provided some more context. Please let us know if you have further questions or suggestions.

We thank the reviewer for the careful and thorough review and provide the explanation below:

• $\bar{C}$ and $\bar{B}$: The notation $\bar{C}$ and $\bar{B}$ above Equation (24) are defined in line 202 and 203, which are our estimates for $C,B$ if the Hankel were exactly known. As the Hankel itself is estimated in the actual Algorithm, Eq. (24) are estimated versions of $\bar{C},\bar{B}$.

• Unknown $k$: What we can prove is that (a) the difference between the lifted Hankel $H$, and $\tilde{H}$ is exponentially small (Lemma C.3); (b) $\tilde{H}$ is formed by a product of a $d_y p\times k$ and a $k\times d_u q$ matrix, so its rank is at most $\min(d_y p, d_u q,k)$; (c) All non-zero singular values of $\tilde{H}$ is exponentially large in $m$. Therefore, when $p,q,m$ are large, from (b) (c) we know $\tilde{H}$ will have $k$ exponentially large (in $m$) singular values and other singular values are zero. Then from (a) we know $H$ will have $k$ exponentially large singular values, and all other singular values are exponentially small. This is the singularvalue gap our paper was referring to and how $k$ can be identified from this gap. For your question where $m,p,q$ are not large enough, there are two cases. Case I: If $p,q$ are not large enough such that the rank of $\tilde{H}$ is smaller than $k$ (i.e. $d_yp<k$ or $d_u q<k$), then by (a) (c) all singular values of $H$ will be exponentially large, and the singular value gap will not be observed. Case II: In the second case, where $m$ is not large enough, then the distinction between the ''exponentially large in $m$'' terms and ``exponentially small in $m$'' will not be obvious, so a large gap is not obvious. To conclude: when $p,q,m$ are not large enough, the singular value gap will not be obvious so one won't mistakenly believe the true $k$ would be identified. Actually, the above discussion also sheds light on how to tune $p,q,m$ to identify the gap.

• Simulation setup: We thank the reviewer for pointing out the issue with the setup in the simulation. In the simulation, $B$ and $C$ are both initialized randomly, where each entry is sampled from an independent uniform distribution between 0 and 1. we still assumed $D = 0$ as discussed in Section 2. Although we discussed how to tackle the case when $D \neq 0$ in Appendix $F$, We did not have non-zero $D$ in the simulation.

• Regarding $\delta$: Yes, if we run the simulation a sufficient number of times for fixed $m,p,q$, we can recover the value of $\delta$. In Theorem 5.3 and Theorem C.2, we use $\delta$ to show that longer trajectories would have better chances of finding accurate Hankel estimation.

We hope the above have answered your questions and provided some more context. Please let us know if you have further questions or suggestions.