Behavioral systems, subspace predictor

1 Introduction↩︎

Prediction from data is a canonical problem across machine learning, signal processing, and control, where a common theme is to exploit low-complexity structure for forecasting and decision making. In systems and control, this viewpoint has led to subspace-based methods for system identification and data-driven prediction/control [1]–[3]. Behavioral systems theory [4], [5] provides a trajectory-based representation of finite-horizon LTI systems: the set of admissible trajectories forms a linear subspace that can be constructed directly from data. The subspace predictor [2] uses this behavioral subspace to predict future outputs from past input-output data and prescribed future inputs.

Robustness and uncertainty quantification for data-driven prediction can be approached either at the data-matrix level [6] or at the behavioral-subspace level [7]. The latter viewpoint has become increasingly relevant in recent work, especially when the object produced by the algorithm is itself a subspace estimate. For example, online subspace estimation methods track time-varying behavioral subspaces [8], [9], while robust least-squares formulations for data-driven control [10] model uncertainty through geometric neighborhoods of a nominal behavior. In these settings, uncertainty is most naturally described by the distance between the true and estimated behavioral subspaces. These developments motivate a coordinate-free sensitivity analysis directly at the level of behaviors. Since subspace uncertainty is naturally measured by distances such as chordal distance, whereas prediction performance is measured by Euclidean output error, it remains necessary to quantify how geometric discrepancies between behaviors affect the resulting predictions.

Recent work [6] derives computable prediction error bounds under data matrix perturbations caused by additive output noise. In contrast, we study perturbations at the level of behavior itself. This geometric viewpoint is more general than an additive-noise model: it can capture uncertainty arising from truncated SVD preprocessing [6], [11], online subspace identification [8], and geometric uncertainty sets used in robust data-driven control [10]. Our goal is to quantify how a geometric error between the true and estimated behavior propagates to the predicted output. When combined with computable chordal distance bounds such as [11] or with subspace-tracking guarantees [8], our result yields prediction-error certificates that are computable from data and noise bounds.

The contribution of this work is twofold. First, we quantify in Theorem 1 the sensitivity of the subspace predictor to behavioral perturbations by deriving an explicit upper bound on the prediction error in terms of the behavioral distance between the ground-truth subspace and its approximation. The bound also depends on a quantitative observability property of the system. Second, we provide a data-dependent prediction-error bound whose constants are computable from the estimated subspace in Theorem 2. Thus, our results convert geometric uncertainty certificates into prediction-error certificates.

The letter is organized as follows: Section 2 formalizes the problem of interest. Section 3 contains the main results (Theorem 1 and Proposition 1). Section 4 empirically studies the sensitivity of the subspace prediction in terms of the behavioral distance through a numerical case study.

Notation: The set of all nonnegative integers is denoted \(\mathbb{Z}_{\geq 0}\). Given \(i, j, T \in \mathbb{Z}_{\geq 0}\) with \(i<j\) and a sequence \(\smash{\{z(t)\}_{t=0}^{T-1}} \subset \mathbb{R}^n\), define \([i,j]=[i, i+1, \dots, j-1, j]\) and \(z_{[i,j]} \mathrel{\vcenter{:}}= [z(i)^\top ~ z(i+1)^\top ~ \cdots ~ z(j)^\top]^\top.\) We use \(\| \cdot \|_2\) to denote both the Euclidean norm for vectors and the spectral norm for matrices, and \(\| \cdot \|_F\) to denote the Frobenius norm. For \(k \leq T\), let \(\mathcal{H}_k(z_{[0,T-1]})\) denote the block-Hankel matrix whose \((\ell+1)^{\text{th}}\) column is \(z_{[\ell,\ell+k-1]}\), for \(\ell=0,\ldots,T-k\). Given \(T \in \mathbb{Z}_{\geq 0}\), the sequence \(z_{[0,T-1]}\) is called persistently exciting of order \(k\) if \(\mathcal{H}_k(z_{[0,T-1]})\) has full row rank. The Moore-Penrose inverse of a matrix \(M\) is denoted by \(M^\dagger\). We use boldface \(\mathbf{U}\) to denote a linear subspace of \(\mathbb{R}^q\), and plain \(U\) to denote a spanning matrix. These notations are connected by \(\mathrm{im~}U = \mathbf{U}\), and will be used interchangeably. The set of orthogonal matrices is denoted \(\mathrm{O}(d)=\{U \in \mathbb{R}^{d \times d}: U^\top U = I\}\). The set of all singular values of a matrix \(U\) is denoted by \(\sigma(U)\). The \(i^\text{th}\) singular value of a matrix \(U\) is denoted by \(\sigma_i(U)\). The diagonal matrix with diagonal entries \(\gamma_1,\cdots,\gamma_d\) is denoted \(\mathop{\mathrm{diag}}(\gamma_1,\cdots,\gamma_d)\). The Grassmannian is defined as \(\mathrm{Gr}(q, d)\mathrel{\vcenter{:}}= \{\mathbf{U}\) is a linear subspace of \(\mathbb{R}^q: \dim \mathbf{U}= d \}.\) Let \(\mathbf{U}, \mathbf{V}\in \mathrm{Gr}(q, d)\) and let \(U, V \in \mathbb{R}^{q \times d}\) be orthonormal bases spanning \(\mathbf{U}\) and \(\mathbf{V}\), respectively. Consider the compact SVD \(U^\top V = P \Sigma Q^\top\), where \(\Sigma=\mathop{\mathrm{diag}}(\sigma_1,\dots,\sigma_d)\) with \(\sigma_1 \geq \dots \geq \sigma_d \geq 0\). The principal angles \(\theta_1, \dots, \theta_d \in [0, \pi/2]\) between \(\mathbf{U}\) and \(\mathbf{V}\) are given by [12] \(\cos \theta_i = \sigma_{i}\). The chordal distance between \(\mathbf{U}, \mathbf{V}\) is defined equivalently by \[\label{eq:chord} d(\mathbf{U}, \mathbf{V}) = \left(\sum_{i=1}^d \sin^2 \theta_i\right)^{\frac{1}{2}} = \frac{1}{\sqrt{2}}\norm{UU^\top - VV^\top}_F.\tag{1}\] Note that the chordal distance is coordinate-free: \(d(\mathrm{im~}(UQ), \mathrm{im~}(VQ)) = d(\mathrm{im~}U, \mathrm{im~}V)\) for all \(Q \in \mathrm{O}(d)\). The Gaussian distribution with mean \(0\) and covariance \(\Sigma \in \mathbb{R}^{d \times d}\) is denoted \(\mathcal{N}(0, \Sigma)\).

2 Problem Statement↩︎

Consider the observable LTI dynamical system \[\label{eq:LTI-1} \begin{align} x(t+1) &= A x(t) + B u(t), \\ y(t) &= C x(t) + D u(t), \end{align}\tag{2}\] where \(A \in\mathbb{R}^{n\times n}\), \(B\in\mathbb{R}^{n\times m}\), \(C\in\mathbb{R}^{p\times n}\), \(D\in\mathbb{R}^{p\times m}\), \(x(t)\) is the state, \(u(t)\) is the input, \(y(t)\) is the output at time \(t \in \mathbb{Z}_{\geq 0}\). The restricted behavior of 2 is defined as \[\begin{align} &\mathcal{B}_{[0,L-1]} = \bigg\{ \left(u_{[0,L-1]}, y_{[0,L-1]} \right) \in \mathbb{R}^{(m+p)L} \bigg|~ \text{there exists} \\ &x_{[0,L]} \in \mathbb{R}^{n(L+1)} \text{ such that \eqref{eq:LTI-1} holds for all t \in [0,L-1]} \bigg\}. \end{align}\] Any \((u_{[0,L-1]}, y_{[0,L-1]}) \in \mathcal{B}_{[0,L-1]}\) is called a trajectory of system 2 . The following lemma [5] shows that \(\mathcal{B}_{[0,L-1]}\) is a linear subspace of \(\mathbb{R}^{(m+p)L}\) and provides a data matrix representation.

Lemma 1 provides the foundation for viewing systems over finite horizons as subspaces [14]. The data-driven prediction problem posed in [2] aims to find the output sequence of an unknown dynamical system corresponding to an input sequence based on a past input-output trajectory collected from the system, without explicitly identifying a state-space model. We now introduce the subspace predictor, which is the key ingredient of data-driven prediction and the main object studied in this work.

Throughout this paper, the subscripts \({u_{\mathrm{p}}}, {u_{\mathrm{f}}}, {y_{\mathrm{p}}}, {y_{\mathrm{f}}}\) always refer to this block partition for any matrices with \((m+p)(T_\mathrm{ini}+T_\mathrm{f})\) rows. When the vector \((u_\mathrm{ini}, u, y_\mathrm{ini})\) is clear from the context, we equivalently write \(y^\mathrm{pred}_{[0,T_\mathrm{f}-1]}(X) \mathrel{\vcenter{:}}= (y^\mathrm{pred}_0(X), \dots, y^\mathrm{pred}_{T_\mathrm{f}-1}(X))\) for \(\mathcal{S}(X, u_\mathrm{ini}, u, y_\mathrm{ini})\). The subspace predictor can be used to predict future output trajectories of a dynamical system given future input \(u \in \mathbb{R}^{mT_\mathrm{f}}\) and an initial trajectory \((u_\mathrm{ini}, y_\mathrm{ini}) \in \mathcal{B}_{[0,T_\mathrm{ini}-1]}\) [2]. We wish to study how the perturbation to \(\mathcal{B}_{[0,T_\mathrm{ini}+T_\mathrm{f}-1]}\) affects the subspace prediction. Let \(\widehat{\mathcal{B}}_{[0,T_\mathrm{ini}+T_\mathrm{f}-1]}\) be the behavior of an approximate system \[\label{eq:LTI-2} \begin{align} \widehat{x}(t+1) &= \widehat{A} \widehat{x}(t) + \widehat{B} u(t), \\ \widehat{y}(t) &= \widehat{C} \widehat{x}(t) + \widehat{D} u(t), \end{align}\tag{3}\] where \(\widehat{A} \in\mathbb{R}^{n\times n}\), \(\widehat{B}\in\mathbb{R}^{n\times m}\), \(\widehat{C}\in\mathbb{R}^{p\times n}\), \(\widehat{D}\in\mathbb{R}^{p\times m}\). Let \(u^\mathrm{d}_{[0,T-1]}\) be persistently exciting of order \(n+T_\mathrm{ini}+T_\mathrm{f}\), let \((u^\mathrm{d}_{[0,T-1]}, y^\mathrm{d}_{[0,T-1]}) \in \mathcal{B}_{[0,T_\mathrm{ini}+T_\mathrm{f}-1]}\) and \((u^\mathrm{d}_{[0,T-1]}, \widehat{y}^\mathrm{d}_{[0,T-1]}) \in \widehat{\mathcal{B}}_{[0,T_\mathrm{ini}+T_\mathrm{f}-1]}\). We arrange them into Hankel matrices \[\label{eq:hankel-data} {!}{ H \mathrel{\vcenter{:}}= \begin{bmatrix} \mathcal{H}_{T_\mathrm{ini}+T_\mathrm{f}}(u^\mathrm{d}_{[0,T-1]}) \\[3pt] \mathcal{H}_{T_\mathrm{ini}+T_\mathrm{f}}(y^\mathrm{d}_{[0,T-1]}) \end{bmatrix}, \widehat{H} \mathrel{\vcenter{:}}= \begin{bmatrix} \mathcal{H}_{T_\mathrm{ini}+T_\mathrm{f}}(u^\mathrm{d}_{[0,T-1]}) \\[3pt] \mathcal{H}_{T_\mathrm{ini}+T_\mathrm{f}}(\widehat{y}^\mathrm{d}_{[0,T-1]}) \end{bmatrix}. }\tag{4}\] Given \((u_\mathrm{ini}, y_\mathrm{ini}) \in \mathcal{B}_{[0, T_\mathrm{ini}- 1]}\) and an input sequence \(u \in \mathbb{R}^{mT_\mathrm{f}}\), we denote the predicted output sequence associated with these Hankel matrices by \({y}^\mathrm{pred}_{[0, T_\mathrm{f}-1]}(H)\) and \(y^\mathrm{pred}_{[0, T_\mathrm{f}-1]}(\widehat{H})\). Note that \((u_\mathrm{ini}, y_\mathrm{ini})\) need not be a \(T_\mathrm{ini}\)-length trajectory of 3 , since the subspace predictor \(\mathcal{S}\) is defined for arbitrary \((u_\mathrm{ini}, y_\mathrm{ini}) \in \mathbb{R}^{(m+p)T_\mathrm{ini}}\). This corresponds to the practical setting in which the approximate subspace, \(\widehat{\mathcal{B}}_{[0,T_\mathrm{ini}+T_\mathrm{f}-1]}\), is used for prediction. Denote \(b=(u_\mathrm{ini}, u, y_\mathrm{ini})\). Denote \(M=[H_{u_{\mathrm{p}}}^\top H_{u_{\mathrm{f}}}^\top H_{y_{\mathrm{p}}}^\top]^\top\) and similar for \(\widehat{M}\), the prediction error can be first bounded as \[\label{eq:first-error-bound} \begin{align} &\norm{y_{[0, T_\mathrm{f}-1]}^\mathrm{pred}(\widehat{H})-y^\mathrm{pred}_{[0, T_\mathrm{f}-1]}(H)}_2 \\ \leq &\left(\norm{\widehat{H}_{y_{\mathrm{f}}}}_2 \norm{\widehat{M}^\dagger - M^\dagger}_2 + \norm{\widehat{H}_{y_{\mathrm{f}}}- H_{y_{\mathrm{f}}}}_F \norm{M^\dagger}_2\right) \norm{b}_2 \\ \leq &\left(\norm{\widehat{H}_{y_{\mathrm{f}}}}_2 \norm{\widehat{M}^\dagger - M^\dagger}_2 + \norm{\widehat{H} - H}_F \norm{M^\dagger}_2\right) \norm{b}_2 \end{align},\tag{5}\] where the last inequality uses the fact that \(\widehat{H}_{y_{\mathrm{f}}}\) and \(H_{y_{\mathrm{f}}}\) are submatrices of \(\widehat{H}\) and \(H\), respectively. However, the bound in 5 is still expressed in terms of matrix norms. As such, it is not intrinsic to the underlying restricted behaviors, and the bound may vary under different matrix representations of the same subspace. Since our goal is to relate prediction error to a geometric discrepancy between restricted behaviors, it is natural to seek a coordinate-free bound formulated directly in terms of \(d( \mathrm{im~}\widehat{H}, \mathrm{im~}H)\), see Fig. 1.

Because each restricted behavior is a linear subspace by Lemma 1, we use chordal distance \(d\) to measure discrepancy between restricted behaviors, and will refer to chordal distance and behavioral distance interchangeably. This gives rise to the central question of this letter:

We address this problem in Section 3. We first examine how observability properties of 2 affect the prediction error bound given by ?? , then we show that the subspace predictor \(\mathcal{S}\) is coordinate-free (Proposition 1). Finally, we relate the Frobenius norm to behavioral distance through a suitable basis representation (Lemma 3), which leads to a novel prediction error bound (Theorem 1).

3 Prediction Error Bound↩︎

We introduce a state-space formulation of the restricted behavior. For \(k \in \mathbb{Z}_{\geq 0}\), define the extended observability matrix and block Toeplitz matrix as \[\mathcal{O}_k = \begin{bmatrix} C \\ CA \\ \vdots \\ CA^{k-1} \end{bmatrix}, \mathcal{T}_k = \begin{bmatrix} D & 0 & \cdots & 0 \\ CB & D & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ CA^{k-2}B & CA^{k-3}B & \cdots & D \end{bmatrix}.\] Let \((u_{[0,L-1]},y_{[0,L-1]}) \in \mathcal{B}_{[0,L-1]}\) with corresponding state sequence \(x_{[0,L]}\). Denote \[\Phi_L \mathrel{\vcenter{:}}= \begin{bmatrix} 0 & I_{mL} \\ \mathcal{O}_L & \mathcal{T}_L \end{bmatrix},\] then \[\label{eq:state-space-behavior} \begin{bmatrix} \mathcal{H}_L\!\left(u_{[0,T-1]}\right) \\[2pt] \mathcal{H}_L\!\left(y_{[0,T-1]}\right) \end{bmatrix} = \Phi_L \begin{bmatrix} \mathcal{H}_1\!\left(x_{[0,T-L]}\right) \\[2pt] \mathcal{H}_L\!\left(u_{[0,T-1]}\right) \end{bmatrix}.\tag{6}\] It is shown in [13] that \(\mathrm{im~}\Phi_L = \mathcal{B}_{[0,L-1]}\). We refer to \(\Phi_L\) as the trajectory generation matrix of 2 . The main theorem below assumes a lower bound on \(\sigma_{\min}(\Phi_L)\), a quantity tied to the observability property of 2 , which is analyzed in Section 3.2. We now state the main theorem.

Theorem 1 answers Problem 1 by quantifying how perturbations of restricted behavior, measured in behavioral distance, affect the prediction error. In particular, the bound in ?? is coordinate-free: it depends only on the subspaces \(\widehat{\mathbf{U}}\) and \(\mathbf{U}\), and not on the particular matrix representations. This shows that sensitivity of subspace prediction can be characterized at the geometric level. The equal-dimensional assumption that \(\dim \widehat{\mathbf{U}} = \dim \mathbf{U}\) in Theorem 1 ensures that the chordal distance is well-defined. This setting is natural for subspace tracking methods with a fixed target dimension [8]. It also covers truncated SVD-type preprocessing [6] when the noisy data matrix is truncated to the prescribed rank. In contrast, if a rank-selection strategy returns a dimension different from the true dimension, then the present chordal-distance bound does not directly apply. Extensions to different-dimensional behavioral subspaces [14] can be handled using distances on the disjoint union of Grassmannians, such as the Grassmannian metric in [15]; we leave this direction for future work. The proof of Theorem 1, presented in Section 3.4, relies on three auxiliary results, developed in the following sections.

3.1 Coordinate-free subspace predictor↩︎

We establish that the subspace predictor depends only on the restricted behavior, rather than on the particular matrix representation.

Since \(M\) has full column rank and \(M\) is a submatrix of \(X\), we have \(\rank X = \rank Y = r\). Consider the QR decompositions [12] \(X=UR_X, Y=VR_Y\), where \(U,V \in \mathbb{R}^{(m+p)(T_\mathrm{ini}+T_\mathrm{f}) \times r}\) are orthonormal and \(R_X, R_Y \in \mathbb{R}^{r \times r}\) are invertible. Since \(\mathrm{im~}X = \mathrm{im~}Y\), we have \(\mathrm{im~}U = \mathrm{im~}V\), hence \(V=UQ\) for some \(Q \in \mathrm{O}(r)\). Denote \(M_U = [U_{u_{\mathrm{p}}}^\top U_{u_{\mathrm{f}}}^\top U_{y_{\mathrm{p}}}^\top]^\top\), then \(M = M_U R_X\), which implies \(M_U\) also has full column rank. Then \((M_U R_X)^\dagger = R_X^{-1} M_U^\dagger\). Therefore, \(y_{[0,T_\mathrm{f}-1]}^\mathrm{pred}(X) = (U_{y_{\mathrm{f}}}R_X)(M_U R_X)^\dagger b = U_{y_{\mathrm{f}}}R_X (R_X^{-1} M_U^\dagger) b = U_{y_{\mathrm{f}}}M_U^\dagger b.\) Similarly, we have \(y_{[0, T_\mathrm{f}-1]}^\mathrm{pred}(Y) = U_{y_{\mathrm{f}}}QR_Y (M_U QR_Y)^\dagger b = U_{y_{\mathrm{f}}}QR_Y R_Y^{-1}(M_U Q)^\dagger b = U_{y_{\mathrm{f}}}QQ^\top M_U^\dagger b = U_{y_{\mathrm{f}}}M_U^\dagger b.\) Thus \(y_{[0, T_\mathrm{f}-1]}^\mathrm{pred}(X) = y_{[0, T_\mathrm{f}-1]}^\mathrm{pred}(Y)\).

Proposition 1 implies that \(\mathcal{S}\) can be viewed as a well-defined mapping on the Grassmannian, i.e., \[\begin{align} \mathcal{S}: \mathrm{Gr}((m+p)(T_\mathrm{ini}+T_\mathrm{f}), r) &\times \mathbb{R}^{m(T_\mathrm{ini}+T_\mathrm{f})+pT_\mathrm{ini}} \to \mathbb{R}^{pT_\mathrm{f}}, \\ \mathcal{S}(\mathbf{U}, u_\mathrm{ini}, u, y_\mathrm{ini}) &= U_{y_{\mathrm{f}}} \begin{bmatrix} U_{u_{\mathrm{p}}}\\ U_{u_{\mathrm{f}}}\\ U_{y_{\mathrm{p}}} \end{bmatrix}^\dagger \begin{bmatrix} u_\mathrm{ini}\\ u \\ y_\mathrm{ini} \end{bmatrix}, \end{align}\] where \(U\) is an orthonormal basis of \(\mathbf{U}\). When the vector \((u_\mathrm{ini}, u, y_\mathrm{ini})\) is clear from context, the notations \[\mathcal{S}(\mathbf{U}), \quad \mathcal{S}(U), \quad y_{[0, T_\mathrm{f}-1]}^\mathrm{pred}(U)\] will be used interchangeably. Now we can replace \(H\) and \(\widehat{H}\) in 5 with their orthonormal bases \(U\) and \(\widehat{U}\), respectively. Then 5 reveals two sources of prediction error: one through the sensitivity of \(M^\dagger\), governed by the singular values of \(M\), and the other is related to the matrix norm \(\|\widehat{U} - U\|_F\) which depends on the basis representations of the restricted behaviors \(\mathrm{im~}\widehat{H}\) and \(\mathrm{im~}H\). We next study the quantitative observability of 2 , which is related to the magnitude of the singular values of \(M\).

3.2 On the quantitative observability condition↩︎

We now turn to the assumption on \(\sigma_{\min}(\Phi_{T_\mathrm{ini}})\) in Theorem 1. The following proposition shows that this assumption is equivalent to observability of system 2 .

Positivity of \(\sigma_{\min}(\Phi_{T_\mathrm{ini}})\) implies that \(\Phi_{T_\mathrm{ini}}\) has full column rank \(n+mT_\mathrm{ini}\). The upper right block \(I_{mT_\mathrm{ini}}\) has rank \(mT_\mathrm{ini}\), so the remaining \(n\) columns must be linearly independent, which is equivalent to \(\rank(\mathcal{O}_{T_\mathrm{ini}})=n\). Conversely, if \(\rank (\mathcal{O}_{T_\mathrm{ini}}) = n\), then the structure of \(\Phi_{T_\mathrm{ini}}\) implies that \(\Phi_{T_\mathrm{ini}}\) has full column rank, and hence \(\sigma_{\min}(\Phi_{T_\mathrm{ini}}) > 0\). This viewpoint also aligns with the quantitative observability condition [16] introduced in the robust fundamental lemma, which assumes a lower bound on \(\sigma_{\min}(\Phi_{T_\mathrm{ini}+T_\mathrm{f}})\). The parameter \(\beta\) quantifies the degree of observability.

Let \(U\) be the orthonormal matrix from the QR factorization \(\Phi_{T_\mathrm{ini}+T_\mathrm{f}} = UR\), where \(R\) is invertible. Then there exists \(Q \in \mathrm{O}(r)\) such that \(U = VQ\). Let \(\Pi\) denote the operator that selects the first \(mT_\mathrm{ini}+mT_\mathrm{f}+pT_\mathrm{ini}\) rows of a matrix such that \(M = \Pi V = \Pi UQ\). Then \(\Pi \Phi_{T_\mathrm{ini}+T_\mathrm{f}} = (\Pi VQ) R = M QR.\) The matrix \(\Pi \Phi_{T_\mathrm{ini}+T_\mathrm{f}}\) admits the form \[\Pi \Phi_{T_\mathrm{ini}+T_\mathrm{f}} = \begin{bmatrix} 0 & I_{mT_{\mathrm{ini}}} & 0\\ 0 & 0 & I_{mT_f}\\ \mathcal{O}_{T_{\mathrm{ini}}} & \mathcal{T}_{T_{\mathrm{ini}}} & 0 \end{bmatrix} \mathrel{\vcenter{:}}= \begin{bmatrix} F ~|~ G \end{bmatrix}\] where \(F=\begin{bmatrix} 0 & I_{mT_\mathrm{ini}} \\ 0 & 0 \\ \mathcal{O}_{T_\mathrm{ini}} & \mathcal{T}_{T_\mathrm{ini}} \end{bmatrix}\) and \(G=\begin{bmatrix} 0 \\ I_{mT_\mathrm{f}} \\ 0 \end{bmatrix}\). Since \(F^\top G = 0\), the column spaces of \(F\) and \(G\) are orthogonal. Therefore for any vector \(x=(x_1, x_2)\) where \(x_1 \in \mathbb{R}^{mT_\mathrm{ini}+pT_\mathrm{ini}}, x_2 \in \mathbb{R}^{mT_\mathrm{f}}\), \[\begin{align} &\| \Pi \Phi_{T_\mathrm{ini}+T_\mathrm{f}} x\|_2^2 = \| F x_1 \|_2^2 + \| G x_2 \|_2^2 \\ \geq &\sigma_{\min}(F)^2\|x_1\|_2^2 + \sigma_{\min}(G)^2\|x_2\|_2^2 \\ \geq &\min(\sigma_{\min}(F)^2, \sigma_{\min}(G)^2) (\|x_1\|_2^2 + \|x_2\|_2^2) \\ = &\min(\sigma_{\min}(F)^2, \sigma_{\min}(G)^2) \|x\|_2^2. \end{align}\] Now \(\sigma_{\min}(G) = 1\), and \(\sigma_{\min}(F) = \sigma_{\min}(\Phi_{T_\mathrm{ini}}) \geq \beta\). Therefore, \(\sigma_{\min}(\Pi \Phi_{T_\mathrm{ini}+T_\mathrm{f}}) \geq \min \{\beta, 1\}.\)

We have the inequality \(\|XYx\|_2 \geq \sigma_{\min}(X)\|Yx\|_2 \geq \sigma_{\min}(X)\sigma_{\min}(Y)\|x\|_2\) for all \(x \in \mathbb{R}^r\), where \(X=\Pi \Phi_{T_\mathrm{ini}+T_\mathrm{f}}, Y=(QR)^{-1}\). Hence \(\sigma_{\min}(XY) \geq \sigma_{\min}(X) \sigma_{\min}(Y)\). Since \(\Phi_{T_\mathrm{ini}+T_\mathrm{f}}\) has full column rank and \(R\) is invertible, \[\begin{align} &\sigma_{\min}(M) = \sigma_{\min}(XY) \geq \sigma_{\min}(X)\sigma_{\min}((QR)^{-1}) \\ &= \sigma_{\min}(\Pi \Phi_{T_\mathrm{ini}+T_\mathrm{f}}) / \sigma_{\max}(QR) \\ &= \sigma_{\min}(\Pi \Phi_{T_\mathrm{ini}+T_\mathrm{f}}) / \sigma_{\max}(R) \geq \min\{1,\beta\} / \alpha. \end{align}\] Finally, since \(\sigma(MP)=\sigma(M)\) for all \(P \in \mathrm{O}(r)\), we have \(\sigma_{\min}(MP) = \sigma_{\min}(M)\). Denote \(\gamma= \min(1, \beta) / \alpha\), then we see that the bound in Theorem 1 improves as \(\beta\) increases, and deteriorates as \(\alpha\) increases. In this sense, a system with stronger quantitative observability leads to a tighter bound, while a larger gain \(\alpha\) leads to greater sensitivity to behavioral perturbations.

3.3 Basis alignment and chordal distance↩︎

The term \(\|\widehat{U} - U\|_F\) in 5 depends on the particular choice of basis of \(\widehat{\mathbf{U}}\) and \(\mathbf{U}\), whereas the prediction error \(\norm{y_{[0,T_\mathrm{f}-1]}^\mathrm{pred}(\widehat{U}) - y_{[0,T_\mathrm{f}-1]}^\mathrm{pred}(U)}_2\) remains invariant by Proposition 1 . It is therefore natural to choose orthonormal bases that minimize \(\|\widehat{U} - U\|_F\). The following lemma is classical and closely related to the orthogonal Procrustes problem [17]. We include it here for completeness, since the proof highlights the geometric role of principal angles and aligns with the subspace-based viewpoint adopted in this work.

Let \(R \in \mathrm{O}(r)\), we have \(\|U - \widehat{U}R\|_F^2 = 2r - 2\operatorname{Tr} (U^\top \widehat{U} R)\). Denoting \(S=Q^\top RP\), we have \(\require{physics} \Tr(U^\top \widehat{U} R) = \sum_{i=1}^r S_{ii} \cos \theta_i.\) Since \(S = Q^\top RP \in \mathrm{O}(r)\), we have \(|S_{ii}| \leq 1\) for \(i=1,\dots, r\). Therefore, \(\operatorname{Tr}(U^\top \widehat{U} R)\) is maximized when \(R = QP^\top\), in which case \(S=I_r\). Then \[{!}{\displaystyle \norm{U - \widehat{U}QP^\top}_F^2 = 2r - 2\sum_{i=1}^r \cos \theta_i \leq 2\sum_{i=1}^r \sin^2 \theta_i \\ = 2d(\widehat{\mathbf{U}}, \mathbf{U})^2, }\] where the last inequality follows from \(\sin^2 \theta_i - (1 - \cos \theta_i) = \cos \theta_i (1 - \cos \theta_i) \geq 0\).

3.4 Proof of Theorem 1↩︎

By Lemma 3, once an orthonormal basis \(U\) of \(\mathbf{U}\) is fixed, we may choose an orthonormal basis \(\widehat{U}\) of \(\widehat{\mathbf{U}}\) so that \(\|\widehat{U}-U\|_F\) is minimized. In the following, \(\widehat{U}\) denotes such a minimizing basis constructed in Lemma 3. Denote \(\kappa = d(\widehat{\mathbf{U}}, \mathbf{U})\), we have \(\|\widehat{U}_{y_{\mathrm{f}}}- U_{y_{\mathrm{f}}}\|_F \leq \|\widehat{U} - U\|_F \leq \sqrt{2}\kappa\). The difference between \(\widehat{M}^\dagger\) and \(M^\dagger\) can be bounded using results on perturbation of pseudo-inverse [18], which implies \[{!}{\displaystyle \norm{\widehat{M}^\dagger - M^\dagger}_2 \leq \frac{1+\sqrt{5}}{2} \max\left\{\norm{\widehat{M}^\dagger}_2^2, \norm{{M}^\dagger}_2^2 \right\} \norm{\widehat{M} - M}_2. }\] By Weyl’s inequality [19] for singular values and Lemma 3, \[{!}{\displaystyle |\sigma_{\min}(\widehat{M}) - \sigma_{\min}(M)| \leq \norm{\widehat{M} - M}_2 \leq \norm{\widehat{M} - M}_F \leq \sqrt{2}\kappa, }\] hence \(\sigma_{\min}(\widehat{M}) \geq \sigma_{\min}({M}) - \sqrt{2}\kappa \geq \gamma- \sqrt{2}\kappa\) by Lemma 2. Since \(\kappa \leq \gamma/ {2\sqrt{2}}\), then \(\sigma_{\min}(\widehat{M}) \geq \gamma/ {2}\), so \(\|\widehat{M}^\dagger\|_2 = \sigma_{\min}(\widehat{M})^{-1} \leq {2} / \gamma\). We then have \[\max\left\{\norm{\widehat{M}^\dagger}_2^2, \norm{{M}^\dagger}_2^2 \right\} \leq \max\left\{\frac{4}{\gamma^2}, \frac{1}{\gamma^2}\right\} = \frac{4}{\gamma^2}.\] Plugging the above inequalities into 5 , using Lemma 3, and using the fact that \(\|\widehat{U}_{y_{\mathrm{f}}}\|_2 \leq \|\widehat{U}\|_2 = 1\), we have \[\begin{align} &\norm{\mathcal{S}(\widehat{\mathbf{U}}) - \mathcal{S}(\mathbf{U})}_2 = \norm{\mathcal{S}(\widehat{U}) - \mathcal{S}(U)}_2 \notag\\ &\leq \left( \frac{2(1+\sqrt{5})\|\widehat{U}_{y_{\mathrm{f}}}\|_2}{\gamma^2} + \frac{1}{\gamma} \right) \norm{\widehat{U}-U}_F\|b\|_2 \\ &\leq \left( \frac{2(1+\sqrt{5})}{\gamma^2} + \frac{1}{\gamma} \right) \sqrt{2}d(\widehat{\mathbf{U}},\mathbf{U})\|b\|_2 . \end{align}\]

The proof uses classical pseudoinverse perturbation tools, but Theorem 1 differs from classical subspace-angle perturbation bounds [20]: it bounds the prediction error of the block-structured subspace predictor under chordal-distance perturbations, with constants containing system-theoretic quantities. Thus, the theorem is a Lipschitz-continuity result for the subspace predictor.

3.5 Data-dependent prediction error bound↩︎

In practice, the prediction error bound ?? is not directly computable because the true behavior is unknown. In this section, we derive a data-dependent error bound, with the only unknown term being the chordal distance between the perturbed and the true subspaces. In some settings, however, this chordal distance can itself be bounded. For noisy data processed by truncated SVD, [11] provides a computable bound for the chordal distance between the estimated and true subspaces. For recursive subspace identification, subspace tracking guarantees such as [8] provide an online bound on the chordal distance of the estimated and true subspaces. Thus, Theorem 2 can be composed with these results to produce computable certificates.

Let \(U\) be an orthonormal basis of \(\mathbf{U}\) and let \(U_{y_{\mathrm{f}}}^{(1)}\) be defined similarly as \(\widehat{U}_{y_{\mathrm{f}}}^{(1)}\). Similar to \(\eqref{eq:first-error-bound}\), we also have \[\begin{align} &\norm{y_0^\mathrm{pred}(\widehat{\mathbf{U}}) - y_0^\mathrm{pred}(\mathbf{U})}_2 \leq \\ &\left( \norm{U^{(1)}_{{y_{\mathrm{f}}}}}_2 \norm{M^\dagger - \widehat{M}^\dagger}_2 + \norm{\widehat{U}^{(1)}_{{y_{\mathrm{f}}}} - U^{(1)}_{{y_{\mathrm{f}}}}}_F \norm{\widehat{M}^\dagger}_2 \right) \|b\|_2. \end{align}\] Let \(\kappa = d(\widehat{\mathbf{U}}, \mathbf{U})\). By Weyl’s inequality used in Section 3.4, we have \(\sigma_{\min}(M)\ge \sigma_{\min}(\widehat M)-\sqrt{2}\kappa \ge {\sigma_{\min}(\widehat M)} / {2},\) so that \(\|M^\dagger\|_2=\sigma_{\min}(M)^{-1} \le {2} / {\sigma_{\min}(\widehat M)}.\) Therefore, \(\max\{\|M^\dagger\|_2^2,\|\widehat{M}^\dagger\|_2^2\} \le {4} / \sigma_{\min} (\widehat{M})^2.\) The rest of the proof proceeds similarly as in Section 3.4, which yields the bound ?? .

4 Numerical Experiment↩︎

We evaluate the data-dependent bound ?? in the Monte Carlo setting of [6]. We generate random stable LTI systems of order \(n\), with \(n\in\{1,2\}\), \(T_\mathrm{ini},T_\mathrm{f}\in\{1,2,3\}\), and \(p,m\) chosen between \(1\) and \(n\). For each random system, we collect \(100\) time steps of offline input-output data using inputs uniformly distributed in \((-1,1)\). The offline output data are then corrupted pointwise by additive noise \(e(t)\) satisfying \(\|e(t)\|_2\leq N\), where \(N\) is the prescribed output-noise bound and \(N \in [10^{-8}, 10^{-3}].\) Further implementation details follow [6]. Code reproducing the experiment is available at github.com/DianJin-Frederick/subspace_prediction_under_behavioral_perturbation.

Our first goal is to compare ?? with the prediction-error bound in [6]. Since the chordal distance \(d(\widehat{\mathbf{U}},\mathbf{U})\) may be unknown in practice, our second goal is to evaluate a fully computable version of ?? by replacing \(d(\widehat{\mathbf{U}}, \mathbf{U})\) with the computable chordal distance bound \[\label{eq:alsalti-bound} d(\widehat{\mathbf{U}}, \mathbf{U}) \leq \frac{C_{\theta}}{\delta_1 \delta_2} N,\tag{7}\] where \(C_\theta, \delta_1, \delta_2\) can be computed from data and are defined in [11]. Substituting 7 into ?? removes the only non-computable term \(d(\widehat{\mathbf{U}},\mathbf{U})\) and yields a prediction-error certificate computable from noisy data and the noise bound \(N\). We use the relative gap defined in [6], which measures the percentage difference between the right-hand side and the left-hand side of the corresponding prediction-error bound, normalized by the output magnitude. Denote \(\mathrm{L} = \norm{y_0^\mathrm{pred}(\widehat{\mathbf{U}}) - y_0^\mathrm{pred}(\mathbf{U})}_2\) and \(\mathrm{R} = \left( \frac{2(1+\sqrt{5})\|\widehat{U}^{(1)}_{y_{\mathrm{f}}}\|_2}{\sigma_{\min}(\widehat{M})^2} + \frac{1}{\sigma_{\min}(\widehat{M})} \right) \sqrt{2}d(\widehat{\mathbf{U}}, \mathbf{U}) \|b\|_2\), then \[\text{relative gap} \mathrel{\vcenter{:}}= \frac{\mathrm{R}-\mathrm{L}}{\|y_0^\mathrm{pred}(\mathbf{U})\|_2}.\]

Results. Fig. 2 shows that ?? , when evaluated with the true chordal distance, has a smaller median relative gap than [6] across all tested noise bounds. Thus, when the behavioral distance is available or accurately estimated, our chordal-distance-based bound gives a tighter prediction-error certificate than the data-matrix perturbation bound of [6]. When the true chordal distance is replaced by the computable bound 7 , the resulting certificate remains fully computable but becomes more conservative. This is expected, since the error is bounded in two stages: first from noisy data to a chordal-distance bound, and then from chordal distance to prediction error. The omitted green outliers correspond to cases where the chordal-distance bound 7 is very conservative [11]. Thus, the experiment illustrates both the tightness advantage of ?? over the baseline and the additional conservatism introduced when the bound is made fully computable.

5 Conclusion↩︎

This work shows two fundamental properties of the subspace predictor. First, we showed that it is coordinate-free, in the sense that it depends only on the restricted behavior of a system. Second, we derived an explicit prediction error bound directly in terms of the chordal distance between two restricted behaviors. This provides a quantitative link between geometric uncertainty and prediction accuracy. Future work includes combining the results with computable subspace-tracking bounds and extending the analysis from prediction to closed-loop data-driven control. We hope that the coordinate-free bound developed will serve as a useful building block for more data-driven control methods.

Acknowledgement↩︎

The authors used ChatGPT to assist with language refinement and editing in parts of this manuscript. All generated content was subsequently reviewed and edited by the authors, who take full responsibility for the final manuscript.

References↩︎