FEDERAL UNIVERSITY OF MINAS GERAIS
Department of Electronic Engineering
Technical Report
Robust Bayesian Subspace Identification for Small Data Sets

Model estimates obtained from traditional subspace identification methods may be subject to significant variance. This elevated variance is aggravated in the cases of large models or of a limited sample size. Common solutions to reduce the effect of variance are regularized estimators, shrinkage estimators and Bayesian estimation. In the current work we investigate the latter two solutions, which have not yet been applied to subspace identification. Our experimental results show that our proposed estimators may reduce the estimation risk up to $40\%$ of that of traditional subspace methods.

In the scope of subspace identification, a number of related works explore regularization techniques such as nuclear norm minimization [1, 2, 3, 4, 5]. Under specific hyperparameter choices, regularized estimators can be understood in a Bayesian framework as maximum a posteriori estimators. These regularized estimators carry two important disadvantages in comparison to our approach: i) the choice of hyperparameters is not straightforward and it often does not take advantage of the problem statistics at hand; ii) they are not oriented toward risk minimization.

The usage of regularized, shrinkage and Bayesian estimators in system identification is discussed in [6, 7, 4] and it is greeted as a promising and refreshing trend in this field. We are only aware of one work where a shrinking function is used for subspace identification [8]. Differently from our approach, the shrinking is applied to all elements of the matrix and not only to the singular values. In a somewhat related application, in [9] we find the application of wavelet thresholding as a preprocessor in a subspace identification setting. Apart from the field of system identification, we also find the application of similar methods in the context of subspace identification in hyperspectral imaging [10].

If we regard subspace identification as a particular application of principal component analysis, our work would be understood as an application of robust principal component analysis, which finds many applications in image and video processing and econometrics [11]. Bayesian solutions to this problem are found in [12, 13, 14]. Our alternating least squares method was conceived as an improvement upon the Bayesian robust principal component analysis of [15] and it is closely related to the tensorial regression present in [16, 17, 18, 19].

A different approach to subspace identification with finite samples is found in [20], where the focus in on deciding the number of samples in order to bound the estimation error with high probability.

Consider the linear state space model in innovation form

where $x_{k}\in\mathbb{R}^{n_{x}}$ is the state, $y_{k}\in\mathbb{R}^{n_{o}}$ is the measured variable, $u_{k}\in\mathbb{R}^{n_{i}}$ is the input variable and $e_{k}$ is the innovations process, assumed to be white Gaussian noise. The model parameters $A,B,C,D$ and $K$ are matrices with the appropriate dimensions. At times, it may be useful to express the same model in the predictor form

where $z_{k}=[u_{k}^{T}~{}~{}y_{k}^{T}]^{T}$, $A_{K}=A-KC$, $B_{K}=[B-KD~{}~{}K]$.

Following the formulation in [21], we make use of the extended state space model

and its predictor form

where the available data is arranged in Hankel matrices defined by

and similarly are defined $U_{f}$ and $E_{f}$. Here $f>n$ is the future horizon and $N$ is a function of the size of the available data set. The state sequence is defined as

The past information is collected in $Z_{p}=[U_{p}^{T}~{}Y_{p}^{T}]^{T}$ up to the horizon $p>n$ and arranged as

and likewise for $Y_{p}$.

As a consequence, $\Gamma_{f}$ is the extended observability matrix as defined by

$H_{f}$ and $G_{f}$ are Toeplitz matrices given by

and

where $\Sigma$ is the innovations covariance. The matrix $H_{fp}$ is reminiscent of the matrix of Markov parameters and is given by

with

for $i=1,2$ and $B_{K}^{(1)}=B-KD$ and $B_{K}^{(2)}=K$. These matrices can also be decomposed as products $H_{fp}^{(i)}=\Gamma_{f}L_{p}^{(i)}$ of the extended observability $\Gamma_{f}$ and controllability matrices:

for $i=1,2$.

Several estimation methods were proposed in the literature departing from this formulation and exploiting the fact that the data as structured in (6) should lie in low dimensional subspaces (see [21, 22] for a comprehensive review). A common approach is to solve (6) via least squares and then estimate $\Gamma_{f}$ from a truncated singular value decomposition of $H_{fp}$. The matrices $A$ and $C$ could then be recovered from estimates of $\Gamma_{f}$ and $\Gamma_{f+1}$. The matrices $B$ and $D$ could then be obtained from a new least squares problem or from estimates of $L_{p}$.

From a general perspective, most methods obtain an initial estimate $\hat{H}_{fp}$ from least squares and, given two weight matrices $W_{1}$ and $W_{2}$, obtain a low-dimensional estimate $\hat{\hat{H}}_{fp}$ by truncating the singular value decomposition of ${\hat{H}}_{fp}$ to the largest $r$ components:

In this context, $r$ would be an estimate of the system order $n$. Methods of obtaining such an estimate of $n$ have been treated with difficulty in the literature and usually one resorts to ad hoc solutions or the visual inspection of the singular values.

Traditional subspace identification methods obtain an initial estimate of $H_{fp}$ from (6) applying the least squares method as follows:

Notice, however, that this estimate only approximates the maximum likelihood estimate as the noise term $G_{f}E_{f}$ in (6) is at best only approximately white and as $H_{f}$ is a lower triangular Toeplitz matrix.

At this point it is worthwhile having a discussion on the approximations typically made in the process of solving (6). Given that $E_{f}$ is structured as a Hankel matrix as in (9), its components are not truly independent. Nevertheless, the off diagonal elements of $\text{E}[E_{f}E_{f}^{T}]$ are reasonably sparse, which gives a reasonable approximation of a white matrix. In order the incorporate the Hankel and Toeplitz structure of $E_{f}$ and $H_{f}$ in the computation, one would have to vectorize equation (6):

where the lower case letters denote the vectorized version of the corresponding matrices. In addition, $h_{f}$ and $e_{f}$ lie in lower dimensional subspaces. This can be solved introducing known matrices $B_{T}$ and $B_{W}$ that convert the original lower dimensional data into the vectorized versions of the Toeplitz matrix $H_{f}$ and the Hankel matrix $E_{f}$. This would result in the new equation

However, the least squares problem that arises in (19) is high dimensional, sparse and poorly conditioned numerically. In this regard, one of the virtues of the formulation in (6) is to avoid the numerical difficulties that would arise in a more detailed formulation.

Most subspace identification methods rely on the limit case of large enough data sets, $N\to\infty$. In addition, most approaches are based on least squares or, equivalently, on maximum likelihood estimators of $H_{fp}$. However, it is well known that such estimators perform poorly in regards to the estimation risk under quadratic loss (are not admissible) when the number of parameters is larger than 2 [23, 24], which is classically known as the Stein phenomenon [25]. On the other hand, the number of parameters in $H_{fp}$ may be as large as $n_{o}\cdot(n_{i}+n_{o})\cdot f\cdot p$, which may be considerably large. Even if we consider that $H_{fp}$ is low rank and may be parameterized and estimated in a smaller space, it would still depend on $(n_{i}+2n_{o})\cdot n_{x}$ parameters, which is still large in most applications.

In addition, there seems to be considerable difficulty in the literature to find appropriate methods to estimate the system order $n_{x}$, which is of particular concern given the high sensitivity of the results to this parameter.

Within this context, our goal is to investigate robustified methods to estimate $H_{fp}$. In particular, two groups of methods shall be considered. The first group is comprised of the main shrinkage estimators of singular values available in the literature. The second group consists of Bayesian estimators that can be efficiently computed with a Gibbs sampling scheme.

Given that the estimation errors in the matrices $(A,B,C,D,K)$ are upper bounded by a constant times $\|H_{fp}-\hat{H}_{fp}\|$ as demonstrated in [20], we concentrate our efforts on the estimators of $H_{fp}$ and leave an investigation on the matrices estimators to future work.

Shrinkage estimators are biased estimators that shrink the maximum likelihood estimate towards a pre-assigned value (often to zero). These estimators were first proposed by Stein [25] and it has been demonstrated that this type of estimator may provide smaller estimation risk than the maximum likelihood estimator under quadratic loss when the number of parameters is larger than 2. The rationale behind such estimators is that, although shrinking introduces bias, it reduces variance and therefore it may be beneficial in reducing the overall quadratic error.

In the context of singular value decompositions, approaches based on random matrix theory obtain optimal shrinkage estimators in the asymptotic case of the matrix dimensions approaching infinity. This assumption may be far from adequate in the case of subspace identification with small data sets. Non-asymptotic results are also available based on the minimization of unbiased estimates.

Most approaches consider the problem of estimating some matrix $X\in\mathbb{R}^{i\times j}$ from a noisy measurement

where $\sigma>0$ is known and the noise matrix $W$ is orthogonally invariant, i.e., its distribution is the same as that of $UWV$ for any orthogonal matrices $U$ and $V$. In practice, $Y$ is not a single measurement, but the maximum likelihood estimate obtained from multiple measurements, which gives a sufficient statistic based on all available measurements.

Considering the Frobenius norm squared as the risk loss function, and given that the noise is orthogonally invariant, it suffices to restrict our search to shrinkage estimators of the form:

where $USV^{T}$ is the singular value decomposition of $Y$ and $\eta(\cdot)$ is the shrinkage function.

Bayesian estimators are inherently robust. Under mild conditions [24], they are admissible, i.e., there exists no estimator that improves their risk for all the parameters in the parameter space. Under group invariance of priors and losses (see [24] for definitions), they are also minimax, i.e., they minimize the risk under the worst case parameter value. The shrinkage estimators presented above may often be characterized as Bayes estimators with an empirical prior, i.e., a prior distribution that is constructed from the data itself.

In our approach to subspace identification, we aim to construct a Bayesian method that is computationally simple. With this in mind, we choose priors that lead to simple regularized least squares steps. The priors are obtained empirically from the data. We have experimented with hierarchical Bayes priors as well but did not find a convincing improvement with respect to the simpler empirical priors. Our model was inspired by that in [15] but, differently from this paper, we do not explicitly estimate the orthogonal matrices $U$ and $V$ and the singular value matrix $S$. We found that better results are obtained by estimating $US^{1/2}$ and $S^{1/2}V^{T}$ instead.

For the noise term $G_{f}$, we propose a prior invariant with respect to the group of lower triangular Toeplitz matrices. Hence, the covariance estimator will be equivariant with respect to this group. In this particular case, it is possible to compute posterior samples in a reasonably simple manner.

Assuming $n_{i}=n_{o}=1$ for simplicity, we adopt the following set of independent priors:

where the $\Xi_{(\cdot)}$ matrices are random matrices with the appropriate dimensions and whose elements are independent normally distributed random variables with variance $1$, the $\Lambda_{(\cdot)}$ matrices are fixed parameters to be specified later and $\bar{G}_{f}=G_{f}/G_{f[1,1]}$. It will be shown later that the improper prior given to $G_{f}$ is invariant under the group of multiplication by lower triangular Toeplitz matrices.

Since the full posterior distribution for the estimation problem at hand is too hard to characterize analytically, we estimate its empirical distribution using a Gibbs sampler. In a Gibbs sampler, samples from dependent variables $x$ and $y$ are drawn iteratively from their conditional distributions as in:

The probability distribution of the resulting Markov chain $(x^{(n)},y^{(n)})$ is shown to converge under mild conditions to $p(x,y)$. The following posterior updates follow from our choice of prior:

and

where $X_{p}^{(n)}=L_{p}^{(n)}Z_{p}$, $\gamma^{(n)}=1/(G_{f[1,1]}^{(n-1)}])^{2}$ and $\Sigma_{e}^{(n)}=G_{f}^{(n-1)}(G_{f}^{(n-1)})^{T}$, the $\Xi_{(\cdot)}^{(n)}$ matrices are random matrices whose components are independently drawn from a unit normal distribution. Intuitively, we are performing independent regressions row by row in (43) and independent regressions column by column in (44), in addition to summing the corresponding simulated noise terms.

In order to define the posterior update for $G_{f}$, we first define the matrix $B_{T}$ that maps the last row of $G_{f}$ to $\text{vec}(G_{f})$, i.e., $\text{vec}(G_{f})=B_{T}(G_{f[i,\cdot]})^{T}$. Similarly, we define the matrix $B_{W}$ such that $\text{vec}(E_{f})=B_{W}e_{[1:i+j-1]}$. Let $\chi_{k}$ denote the chi distribution with $k$ degrees of freedom, then the posterior update for $G_{f}$ is given by:

where $\Omega_{L}^{(n)}$ denotes the lower triangular part of the Cholesky decomposition of $\Omega^{(n)}$. In other words, we first obtain a sample of the last row of $G_{f}^{-1}$, then we construct the full lower triangular Toeplitz matrix $G_{f}^{-1}$ from this row and compute its inverse. Given that the matrix $B_{W}$ is very large and sparse, (46) may be somewhat tricky to compute. A safer alternative may come from ignoring the Hankel structure of $E_{f}$ and assuming that its elements are mutually independent. In this case, the posterior update becomes

Finally, the estimate for ${H}_{fp}$ is obtained by averaging the over the trajectory of the Markov chain:

where $N_{o}$ is some burn-in period intended to remove the effect of transients. In order to reduce variance and improve convergence, one may prefer to average over the expected values of (43) and (44) in every step (obtained by setting the respective $\Xi_{(\cdot)}$ matrices to zero):

Regarding the parameters in the priors, we initially obtain estimates ${H}_{fp}^{(1)}$ and ${H}_{f}^{(1)}$ from (17) and next some truncated singular value decomposition $H_{fp}^{(1)}Z_{p}\approx U_{r}S_{r}V_{r}^{T}$. Then, we make $G_{f}^{(1)}=I_{i}$ and

To justify such a choice, we note that the priors (36) and (38) approximately describe an SVD decomposition:

where ${\Xi_{\Gamma}}/{i}$ and ${\Xi_{L}}/{j}$ behave as orthogonal matrices in expectation, i.e., $\operatorname{E}[\Xi_{\Gamma}^{T}\Xi_{\Gamma}/{i^{2}}]=\operatorname{E}[\Xi_{L}\Xi_{L}^{T}/{j^{2}}]=I_{r}$.

In order to test our proposed methods, we ran Monte Carlo simulations on a large number of systems and compared the estimation risks.

The system order $n_{x}$ was uniformly distributed from $1$ to $10$ and we fixed $n_{i}=n_{o}=1$. The sample size was defined as $N=\lfloor 80\sqrt{n_{x}}\rfloor$. The row-length of Hankel matrices was set as $i=\lfloor N/10\rfloor$.

To obtain a stable system matrix $A$, we first sampled an auxiliary matrix $\tilde{A}$ such that $\tilde{A}_{[k,l]}\sim N(0,1)$ and a spectral radius $\lambda_{a}\sim\mathcal{U}(0,1)$. Then we made $A=\tilde{A}/\lambda_{\max}(\tilde{A})\cdot\lambda_{a}$. The matrices $B$ and $C$ were sampled such that $B_{[k,l]},C_{[k,l]}\sim N(0,1)$. We make $D=0$. The measurement noise covariance $R_{v}$ was generated such that $(R_{v})^{1/2}_{[k,l]}\sim N(0,1)$. Likewise, the process noise covariance $R_{w}$ was generated such that $(R_{w})^{1/2}_{[k,l]}\sim N(0,1)$. The Kalman gain $K$ was therefore obtained from the previous parameters.

The system input $u_{k}$ is comprised of independent samples from $N(0,\textsf{SNR})$, where the signal to noise ratio SNR was sampled uniformly on logarithmic scale such that $\log_{10}\textsf{SNR}\sim\mathcal{U}(-1,2)$.

We applied the proposed algorithms to each model realization and response realization. For each realization and each estimator, we computed the risk

Given that the system models have different scales in each realization, we normalized the risk performance $\mathcal{R}_{k}$ by that of a reference estimator $\mathcal{R}_{k}^{o}$ and averaged over N realizations as such:

The logarithmic function is used to give higher and symmetric weights to risks that are either too low or too high compared to the reference.

We used the same estimate of $\hat{G}_{f}$ for all shrinkage methods as defined by (35). We started by truncating $\hat{H}_{fp}$ to the largest $r$ singular values and then, computing the corresponding $\hat{G}_{f}(r)$ and hard threshold $\lambda_{\text{soft}}(r)$, we obtained

In words, $r^{*}$ is the least $r$ such that the order estimate is less than $r$. We then make $\hat{G}_{f}=\hat{G}_{f}(r^{*})$ and use $r^{*}$ in (35).