11institutetext: University of Hildesheim, 31141 Hildesheim, Germany
22institutetext: Volkswagen AG, Berliner Ring 2, 38440 Wolfsburg, Germany
22email<EMAIL_ADDRESS>
# Towards Comparable Knowledge Distillation in Semantic Image Segmentation

Onno Niemann 11 0009-0004-6131-529X Christopher Vox 22 0009-0003-4050-5070
Thorben Werner 11 0000-0002-5944-6045

###### Abstract

Knowledge Distillation (KD) is one proposed solution to large model sizes and
slow inference speed in semantic segmentation. In our research we identify 25
proposed distillation loss terms from 14 publications in the last 4 years.
Unfortunately, a comparison of terms based on published results is often
impossible, because of differences in training configurations. A good
illustration of this problem is the comparison of two publications from 2022.
Using the same models and dataset, Structural and Statistical Texture
Distillation (SSTKD) [14] reports an increase of student mIoU of 4.54 and a
final performance of 29.19, while Adaptive Perspective Distillation (APD) [30]
only improves student performance by 2.06 percentage points, but achieves a
final performance of 39.25. The reason for such extreme differences is often a
suboptimal choice of hyperparameters and a resulting underperformance of the
student model used as reference point.

In our work, we reveal problems of insufficient hyperparameter tuning by
showing that distillation improvements of two widely accepted frameworks, SKD
and IFVD, vanish when hyperparameters are optimized sufficiently. To improve
comparability of future research in the field, we establish a solid baseline
for three datasets and two student models and provide extensive information on
hyperparameter tuning. We find that only two out of eight techniques can
compete with our simple baseline on the ADE20K dataset.

###### Keywords:

Knowledge Distillation Efficient Semantic Segmentation Model Compression.

mIoU

    Mean Intersection over Union

##  1 Introduction

Advances in Deep Learning techniques for Semantic Image Segmentation brought
major performance improvements to fields such as autonomous driving, medical
image analysis, robotic perception and video surveillance. As these
performance gains often came at the price of increased model complexity and
required computational power, efficient deep learning techniques have become
increasingly relevant [18].

Two techniques, which start with a large model and make it more efficient are
model Pruning and Quantization. Pruning shrinks a model by dropping less
important nodes and Quantization reduces the numerical precision of weights.
Knowledge distillation takes a different approach and does not change model
efficiency during training. Instead it starts with a small model (student) and
improves its performance by leveraging guidance from a more complex model
(teacher). Teacher model weights are frozen and knowledge is distilled into
the student by making an addition to the student loss, which penalizes
differences in student and teacher output. Incentivizing the student to mimic
the more complex behaviour of the teacher can significantly lift its
performance.

KD was introduced in the image classification domain and originally the
distillation loss was applied to student and teacher only at output-level. The
extension of teacher guidance to intermediate layers and the transfer to
semantic segmentation have been studied in various publications. Re-using the
most basic, output-level distillation loss in a segmentation context is
straightforward, as it can be applied on pixel instead of image-level in the
student and teacher output. However, several papers criticize that this naive
approach treats pixels in isolation and introduce more complex techniques.
Most of them build upon the naive pixel-wise distillation by adding loss terms
to it. A problem with these published methods is that a significant share
focuses solely on performance improvement in their own training framework and
fails to provide a good baseline for comparison against other literature.

The dimension of this comparability problem is illustrated by the overview
over important publications and their results we provide in Table 1. The
comparison clearly shows that the baseline performance of the student model
varies strongly between publications, making it hard to compare the quality of
the proposed techniques based on the lift they provide. To address this
problem of comparability, we perform an extensive hyperparameter optimization
of the most fundamental of KD frameworks for semantic segmentation and provide
detailed information on optimal hyperparameters for training of two student
models on three datasets. As a byproduct of this optimization, we find that
the temperature parameter used to ”soften” student and teacher output in image
classification can improve distillation in segmentation, although many
publications in the field ignore it.

In summary, we point out challenges in comparing different methods and present
a training procedure, which sets the ground for a fair comparison of KD
techniques. We further put the performance of three commonly used loss terms
into perspective by comparing them to our achieved results.

##  2 Related Work

#### Semantic Segmentation.

Most earlier image segmentation techniques were based on Partial Differential
Equations or Random Forest methods [10], before the advent of deep learning
sparked a series of publications leveraging the powers of Convolutional Neural
Networks (CNNs). The majority of this early CNN-based segmentation research
focusses on improving model performance, which is achieved by adding skip
connections [17] and decoder networks [19, 24, 2], or concatenating
convolutions at different scales in the pyramid pooling module of Pyramid
Scence Parsing Networks (PSPNets) [44]. More recent publications investigate
the application of transformer based models [45, 28, 33]. Another branch of
segmentation research focuses more on model efficiency instead of performance.
Real-time semantic segmentation aims for fast inference speed while
maintaining performance [43, 21, 15, 37], but there is always a trade-off
between performance and efficiency.

#### Knowledge Distillation in Classification.

KD aims to increase the performance of a compact student model by leveraging
guidance from a large teacher model during training. It is one of the most
widely used model compression techniques due to its broad applicability [31].
Unlike Model Pruning or Quantization, KD is model agnostic [3] as no
requirements are imposed on the teacher or student model. Additionally, KD
allows leveraging unlabelled data as the teacher can provide soft labels for
student training [31].

Generally, the teacher model weights are frozen during student training and
the student is encouraged to mimic the teacher’s output logits by training on
a weighted combination of standard cross entropy and KD loss [4, 13] (Eq. 1).
The cross entropy term (Eq. 2) ensures performance on the labeled training
data while the distillation term (Eq. 3) penalizes deviations from the teacher
output. The weight of the distillation loss, $\lambda$, is a hyperparameter of
student training. Scaling output logits of both models by $\tau>1$ is crucial
for a successfull distillation of knowledge in image classification [13] and
extending the level of student teacher matching to intermediate feature layers
can give a further boost [23, 40].

$L_{S}=L_{CE}+\lambda*L_{KD}$ (1)

$L_{CE}(z_{s},y)=-\sum_{c=1}^{C}y_{c}\log\sigma_{c}(z_{s})$ (2)

$L_{KD}(z_{s},z_{t})=-\tau^{2}\sum_{c=1}^{C}\sigma_{c}(\frac{z_{t}}{\tau})\log\sigma_{c}(\frac{z_{s}}{\tau})$
(3)

Despite the simplicity of the classical KD framework, the underlying
mechanisms are still not well-understood. Commonly, the success of KD is
associated with the information contained in teacher predictions for wrong
classes, referred to as ”dark knowledge” [39]. However, phenomena such as
students outperforming teachers when trained solely on teacher output and high
disagreement between student and teacher predictions suggest that there are
other explanations for the success of KD [27]. Two alternative explanations
are that teacher guidance has a regularizing effect similar to label smoothing
[39] and that teacher output provides a sample-wise importance weighting,
making the student focus on samples of low teacher confidence [8].

#### KD in Image Segmentation.

Applying the classical KD framework to semantic image segmentation problems is
straight-forward. The distillation loss in Equation 3 is usually applied on
image-level, summing over all classes, but since image segmentation is
essentially pixel-wise classification, it can be applied to each pixel in the
image instead. This most basic distillation loss is referred to as pixel-wise
distillation ($L_{PI}$) and is a useful baseline for more complex distillation
schemes. Since pixel-wise distillation treats each pixel in isolation and
ignores the fact that segmentation depends strongly on contextual information,
various alternative distillation techniques exist, many of which use pixel-
wise distillation as one part of their framework.

The earliest proposed technique investigates a ”consistency” loss comparing
regional differences in student and teacher output [34] aiming for a more
contextual distillation of knowledge by matching the distance of the center
pixel and an 8-neighborhood. Instead of a direct matching of student and
teacher output, in Knowledge Adaptation (KA) [12] teacher output is compressed
to a more dense latent space by an autoencoder before being compared to
student logits. Additionally, an affinity loss term is introduced to better
capture long-range dependencies.

Structured Knowledge Distillation (SKD) [16] is the most cited publication in
the field and uses a combination of three loss terms, one of them the basic
pixel-wise loss. Again the other two loss terms are supposed to focus more on
contextual information in both intermediate and output layers. The pair-wise
loss ($L_{PA}$) is based on a pair-wise Markov random field framework and
encourages students to mimic the teacher at intermediate layers. The holistic
loss ($L_{HO}$) requires an additional discriminator model with the task to
differentiate student and teacher output. Student and discriminator compete
with each other following a similar training protocol as used to train
Generative Adversarial Networks[9], ultimately leading to the student output
being as similar as possible to the teacher output [16].

Building upon SKD, Intra-Class Feature Variation Distillation (IFVD) [32] also
uses the pixel-wise and holistic loss terms, but replaces the pair-wise loss
with the new IFV loss, $L_{IFV}$. Per-class prototypes are calculated and the
distance of intermediate features to the respective prototype is aligned
between student and teacher. Another suggested approach, CSCACE [20],
introduces a Channel and Spatial Correlation (CSC) and an Adaptive Cross
Entropy (ACE) term. CSC calculates correlation matrices between intermediate
features and can be understood as an extension of the pair-wise loss of SKD.
ACE combines the classic pixel-wise KD loss with the ground truth label by
only using the teacher output when the teacher prediction is correct.

As one of the few approaches that modify the pixel-wise distillation loss,
Channel-Wise Distillation (CWD) [26] proposes a channel/class-wise
normalization of student and teacher outputs before calculating the KL-
divergence, a technique that is copied by several later methods. Double
Similarity Distillation (DSD) [7] introduces a pixel-wise similarity loss,
which matches intermediate student and teacher features by utilizing self-
attention maps across multiple layers. A category-wise similarity distillation
loss makes the student mimic the teacher at output level by minimizing the L2
distance between correlation matrices of student and teacher output. Masked
Generative Distillation (MGD) [36] is a more general technique, which can also
be used for knowledge distillation in image classification or object
detection. The proposed method masks some parts of the input while still
requiring the student to mimic the full teacher output.

Inter-Class Distance Distillation (IDD) [41] and Feature-Augmented Knowledge
Distillation (FAKD) [38] use the pixel-wise-distillation loss with channel-
wise normalization as suggested in CWD [26]. IDD additionally includes all
terms introduced by SKD and an inter-class feature distance loss following a
similar reasoning as IFVD. FAKD does not introduce additional loss terms but
proposes the perturbation of intermediate student features in multiple ways
and training the student to mimic the teacher despite the applied
perturbations.

Structural and Statistical Texture Knowledge Distillation (SSTKD) [14] and
Adaptive Perspective Distillation (APD) [30] use the pixel-wise distillation
loss. SSTKD additionally includes the holistic loss of SKD and two novel loss
terms, which encourage the student to mimic low-level texture information of
the teacher [14]. The authors of APD argue that segmentation networks learn to
generalize and thus acquire a universal perception. Their introduced adaptive
perspective includes calculating class-wise average features of individual
images. According to the authors, this process distills contextual information
more explicitly. APD is the only reviewed technique that updates part of the
teacher model during student training.

Self-attention and Self-distillation [1] introduces a self-attention loss term
to make the student learn contextual information from the teacher and a layer-
wise context distillation loss. The second term is different from other
discussed techniques, as it is applied across student layers to ensure a
consistent representation of contextual information in shallow layers. Cross-
Image Relational Knowledge Distillation (CIRKD) [35] again uses the pixel-wise
distillation loss and uses three more loss terms. Unique about this method is
the introduction of a pixel queue to make student models mimic the teachers
distance to output for pixels of the same class from previous images.

##  3 Related Work Comparison

A comparison of the discussed knowledge distillation frameworks is rarely
straightforward. A very common problem is differences in student and teacher
architectures. Especially earlier methods use a variety of models and results
cannot be compared to recent literature. Table 1 presents the results of
different distillation frameworks that report performance on the most common
choice of student and teacher model, PSPNets with ResNet101 and ResNet18
backbone, respectively. This choice of models has been the standard in the
field since the publication of SKD [16], but Consistency [34], CSCACE [20], KA
[12] and SALC [1] use other architectures and are excluded from our table.

Table 1: Comparison of reported student mIoU of various techniques following the most common architectural setup. The ”Ls” column refers to the number of loss terms in addition to the CE loss. ”S only” shows the performance of the student model trained without teacher. ”+T ”contains the best performance reported in the publication applying all proposed loss terms. Methods followed by * use the $L_{PI}$, the ones followed by ${\dagger}$ additionally do channel-wise normalization before applying $L_{PI}$. |  | ADE20K | Cityscapes
---|---|---|---
| Ls | S only | \+ T | S only | \+ T
SKD[16]* | 3 | 33.82 | 36.55 | 69.10 | 72.67
IFVD[32]* | 3 | - | - | 69.10 | 74.54
CWD[26]${\dagger}$ | 1 | 24.65 | 26.80 | 70.09 | 75.90
DSD[7] | 2 | 33.80 | 38.00 | 69.42 | 73.20
MGD[36] | 1 | - | - | 69.85 | 74.10
IDD[41]${\dagger}$ | 6 | 24.65 | 27.69 | 70.09 | 77.59
FAKD[38]${\dagger}$ | 1 | 29.42 | 35.30 | 68.99 | 74.75
SSTKD[14]* | 4 | 24.65 | 29.19 | 69.10 | 75.15
APD[30]* | 3 | 37.19 | 39.25 | 74.15 | 75.68
CIRKD[35] | 4 | - | - | 72.55 | 74.73

Even though the techniques in Table 1 use the same model architectures we can
observe significant performance differences for both the student-only training
(”S only”) and the full KD framework (”+T”). While the differences in ”+T” are
to be expected for different KD algorithms, the inconsistent baseline
performance of the student model makes it very hard to compare different
algorithms.

Additionally, the listed techniques often propose a number of loss terms
(”Ls”), but do not provide comprehensive ablation studies that isolate the
effect of individual terms across different datasets. Finally, many results in
recent literature are reported without standard deviations, with SKD being the
only exemption in Table 1. This introduces another complicating factor when
comparing results of KD algorithms.

For a specific example we look at APD vs. DSD: Judging only by the final
performance of 39.25, APD outperforms all other techniques on ADE20K.
Considering its starting point of 37.19, however, student performance is only
lifted by 2.06. As DSD achieves a lift of 4.20 using one loss term less, but
only achieves a final performance of 38.00, it is unclear, which of the two
techniques is superior.

##  4 Methodology

In this work, we present a series of experiments highlighting the importance
of individual hyperparameter tuning for every combination of datasets and
models. We grid search initial learning rate $\mu(0)$ and regularization rate
$\gamma$ for student with and without teacher model separately to ensure the
best performance possible in both settings. In this way we guarantee a fair
performance comparison when measuring the lift from KD later. In a second
step, we optimize the temperature parameter $\tau$ of the combined student and
teacher training.

We also fine-tune the weight of the pixel-wise distillation loss,
$\lambda_{PI}$, and find an optimal performance for $\lambda_{PI}=1e-1$ for
all experiments. Finally, we test different proposed loss terms in isolation
to measure their effect on student performance by optimizing their respective
loss weights.

##  5 Experiments

###  5.1 Datasets and Evaluation Metric

We evaluate our approach on three public semantic segmentation benchmark
datasets, PascalVOC[6], Cityscapes[5], and ADE20K[46]. All datasets have
pixel-wise annotations.

Cityscapes is an urban street scene understanding dataset showing street
scenes recorded in 50 cities. It contains 5,000 finely annotated images with
labels from 19 classes. All images have dimension 2048x1024 and train,
validation and test set contain 2,975, 500 and 1,525 images, respectively.

ADE20K is a complex scene understanding dataset containing 20K/2K/3K images of
different sizes in train, validation and test set. Images show objects, parts
of objects and stuff from varying context and pixels are assigned one of 150
object and stuff class labels.

PascalVOC contains 1,464 images for training, 1,449 images for validation and
a private test set. It has 20 different object classes and images are of
varying size.

Mean Intersection over Union is the metric used in all our experiments.

###  5.2 Implementation Details

As most other relevant methods (Table 1), we follow the architectural setup
and training procedure of SKD [16]. To be precise, this means using PSPNets
[44] with different backbones as student and teacher models. The teacher model
has a ResNet101 [11] backbone and for the student we test a ResNet18[11] and
an EfficientNet-B0 [29] backbone. For simplicity we will refer to the two
student models as ResNet and EffNet students. The learning rate $\mu(i)$
decays over training steps according to Eq. 4.

$\mu(i)=\mu(0)*(1-\frac{i}{\eta})^{0.9}\hskip 1.42262pt;\hskip
5.69054pti:=[1\ldots\eta]$ (4)

$\eta$ is the total number of training batches and i is the current batch.
Unless otherwise stated the ResNet and EffNet student backbones are
initialized with weights pre-trained on ImageNet. The pre-trained student
weights were obtained from the torchvision package (v0.12) of the PyTorch
library [22] for both backbones and the teacher weights were taken from [16]
for Cityscapes and from [42] for the other datasets.

All experiments are conducted with a batch size of 8 and students are trained
with crops of size $512\times 512$ on Cityscapes and $473\times 473$ on ADE20K
and PascalVOC.

All experimental results are calculated on the validation datasets of
Cityscapes, PascalVOC and ADE20K.

###  5.3 The Impact of Temperature

The authors of the original KD framework [13] strongly emphasize the
importance of logit scaling by a temperature parameter $\tau>1$ in image
classification, the main reason for its success being the ”softening” of
teacher output class distributions [13]. To show that teacher output
distributions are ”hard” also in image segmentation we analyze the effect of
different values of $\tau$ on the Shannon Entropy [25] of teacher output. When
all probability mass is assigned to one class the Shannon Entropy is zero and
when the probability mass is distributed evenly over all classes Shannon
Entropy is maximal.

We generate teacher output for 800 randomly selected images from the
Cityscapes dataset, resulting in probability distributions over the 19 classes
for 209,715,200 pixels. These distributions are scaled with different
temperature values according to Equation 3. The shares of Shannon Entropies
over all pixels for different temperatures are shown in Figure 1. Temperatures
of 1, 2, 4, 8, and 16 were chosen as [13] report successful distillation for
values up to 20.

Figure 1: Effect of logit scaling on Shannon Entropy of teacher output
probability distributions over classes. The histograms are based on class
probability distributions of randomly selected pixels from Cityscapes. Without
scaling ($\tau=1$) over 60% of pixels are assigned one class with a
probability close to 1.

An important observation is that when the teacher output is not scaled
($\tau=1$), the entropy strongly spikes at 0. More than 60% of pixels are
classified with a confidence close to 1, suggesting that teacher output might
be too ”hard” for efficient distillation of knowledge and higher temperatures
might help. On the other hand, for a value of $\tau=16$ almost 100% of
distributions have an entropy of approximately 2.9, indicating an almost even
distribution of probability mass over all classes.

###  5.4 Hyperparameter Optimization

As described previously, we optimize hyperparameters separately for student
only and the combined student and teacher training. The results of the grid
search for initial learning rate $\mu(0)$ and weight decay $\gamma$ for both
students on Cityscapes are presented in Table 2(a), right of it in Table 2(b)
are the optimal hyperparameters for all datasets.

Table 2: Grid search results for student-only training. a) shows the whole
grid for Cityscapes, b) the best hyperparameters for all datasets. The best
performance in a) is highlighted in bold.

(a) Cityscapes Grid Search | $\mu(0)$
---|---
| EffNet | ResNet
$\gamma$ | 1e-1 | 5e-2 | 1e-2 | 5e-3 | 5e-2 | 1e-2 | 5e-3
5e-4 | 35.0 | 43.9 | 59.8 | 60.1 | 68.89 | 70.55 | 67.95
5e-5 | 62.2 | 63.6 | 61.1 | 60.1 | 69.74 | 69.48 | 67.98
5e-6 | 64.6 | 63.6 | 61.1 | 59.4 | 70.55 | 67.65 | 67.09

(b) All Datasets | EffNet | ResNet
---|---|---
| $\mu(0)$ | $\gamma$ | $\mu(0)$ | $\gamma$
PascalVOC | 1e-2 | 5e-4 | 5e-3 | 5e-4
Cityscapes | 1e-1 | 5e-6 | 1e-2 | 5e-4
ADE20K | 5e-3 | 5e-5 | 1e-2 | 5e-5

We use the same grid for the student and teacher training and set the
temperature parameter to $\tau=1$. The results are shown in Tables 3 where
again the left shows the detailed grid for Cityscapes and the right the
optimal hyperparameters for all students and datasets.

Table 3: Grid search results for student and teacher training. a) shows the
whole grid for Cityscapes, b) the best hyperparameters for all datasets. The
best performance in a) is highlighted in bold.

(a) Cityscapes Grid Search | $\mu(0)$
---|---
| EffNet | ResNet
$\gamma$ | 1e-1 | 5e-2 | 1e-2 | 5e-3 | 5e-2 | 1e-2 | 5e-3
5e-4 | 40.4 | 46.9 | 62.8 | 61.1 | 66.44 | 71.64 | 70.54
5e-5 | 62.1 | 65.6 | 63.2 | 62.0 | 67.19 | 70.82 | 68.42
5e-6 | 65.0 | 66.1 | 63.1 | 61.6 | 70.02 | 69.72 | 68.64

(b) All Datasets | EffNet | ResNet
---|---|---
| $\mu(0)$ | $\gamma$ | $\mu(0)$ | $\gamma$
PascalVOC | 1e-2 | 5e-6 | 5e-3 | 5e-5
Cityscapes | 5e-2 | 5e-6 | 1e-2 | 5e-4
ADE20K | 1e-2 | 5e-5 | 1e-2 | 5e-5

In a second stage of hyperparameter optimization, we tune the temperature
parameter $\tau$. The tested values and their responses for all models and
datasets can be found in Table 4. A choice of $\tau=1$ appears to work well on
PascalVOC and ADE20K, but greater values yield performance gains on
Cityscapes.

Table 4: Impact of the temperature parameter $\tau$ on segmentation performance. The best performances are highlighted in bold. Results within one standard deviation of the best result are underlined. For ADE20K only one run was computed. | PascalVOC | Cityscapes | ADE20K
---|---|---|---
$\tau$ | EffNet | ResNet | EffNet | ResNet | EffNet | ResNet
1 | 66.24 $\pm$ 0.46 | 65.22 $\pm$ 0.44 | 65.26 $\pm$ 0.59 | 71.33 $\pm$ 0.85 | 36.32 | 37.74
2 | 66.08 $\pm$ 0.16 | 64.37 $\pm$ 0.08 | 64.62 $\pm$ 0.77 | 72.15 $\pm$ 0.10 | 35.51 | 36.85
3 | 65.91 $\pm$ 0.26 | 63.30 $\pm$ 0.44 | 65.08 $\pm$ 0.63 |  72.75 $\pm$ 0.29 | 34.87 | 36.04
4 | 65.02 $\pm$ 0.82 | 63.42 $\pm$ 0.54 | 65.26 $\pm$ 0.43 | 72.44 $\pm$ 0.31 | 35.02 | 36.34
6 | 66.11 $\pm$ 0.67 | 63.65 $\pm$ 0.51 | 64.52 $\pm$ 0.45 | 72.49 $\pm$ 0.43 | 35.23 | 36.49
8 | 66.11 $\pm$ 0.13 | 64.18 $\pm$ 0.36 | 65.58 $\pm$ 0.03 | 72.09 $\pm$ 0.50 | 34.91 | 36.51

###  5.5 Additional Losses

The Losses $L_{PA}$ [16], $L_{HO}$ [16], and $L_{IFV}$ [32] are part of many
KD algorithms, but their individual contribution to the student’s learning
behavior when hyperparameters are optimal is unclear in literature. We use the
optimal hyperparameters of the combined training in Table 3(b) and optimize
the weight of each individual loss term $\\{\
\lambda_{PA},\lambda_{HO},\lambda_{IFV}\\}$ to maximize performance.

Table 5: Impact of individual loss terms on segmentation performance on
Cityscapes. The results for three loss terms, $L_{PA}$, $L_{HO}$ and $L_{IFV}$
are presented in tables a), b) and c), respectively. The best mIoU is
highlighted in bold for each student model.

(a) $L_{PA}$ | EffNet | ResNet
---|---|---
$\lambda_{PA}$ |  |
1e-3 | 65.08 | 70.51
1e-2 | 65.17 | 70.28
1e-1 | 63.09 | 72.98
1e+0 | 13.49 | 71.88
1e+1 | 3.211 | 69.34

(b) $L_{HO}$ | EffNet | ResNet
---|---|---
$\lambda_{HO}$ |  |
1e-4 | 65.84 | 70.35
1e-3 | 65.48 | 70.81
1e-2 | 64.70 | 72.50
1e-1 | 41.53 | 28.20

(c) $L_{IFV}$ | EffNet | ResNet
---|---|---
$\lambda_{IFV}$ |  |
1e-4 | 65.78 | 70.25
1e-3 | 66.01 | 71.70
1e-2 | 65.80 | 71.75
1e-1 | 64.31 | 71.48
1e+0 | 65.03 | 71.75
1e+1 | 65.52 | 70.60
5e+1 | 64.76 | 71.21
1e+2 | 62.70 | 70.20

The results for the Cityscapes dataset are presented in Table 5.

For the EffNet student the optimal values of loss weights are consistently
smaller than for the ResNet student.

###  5.6 Final Performance Comparison

Table 6 shows the results of the additional loss experiments in comparison to
student-only or simple pixel-wise distillation training. It is clear that
adding the teacher model and the pixel-wise distillation loss, $L_{PI}$,
improves student performance, while none of the three tested additional loss
terms provides a further lift. An exception is the PascalVOC dataset, where
the conclusion is less clear.

Table 6: Evaluation of more complex loss terms $L_{PA}$, $L_{HO}$, and $L_{IFV}$ when added to $L_{CE}+L_{PI}$ with optimal hyperparameters. Best results are highlighted in bold; results within one standard deviation of the best result are underlined. | Pascal Voc | Cityscapes | ADE20K
---|---|---|---
$L_{PI}$ | $L_{PA}$ | $L_{HO}$ | $L_{IFV}$ | EffNet | ResNet | EffNet | ResNet | EffNet | ResNet
|  |  |  | 65.43 $\pm$ 0.38 | 64.50 $\pm$ 0.34 | 64.42 $\pm$ 0.46 | 70.45 $\pm$ 0.10 | 34.10 | 35.23
x |  |  |  | 66.24 $\pm$ 0.46 | 65.22 $\pm$ 0.44 | 65.58 $\pm$ 0.03 | 72.75 $\pm$ 0.29 | 36.32 | 37.74
x | x |  |  | 66.40 $\pm$ 0.35 | 65.84 $\pm$ 0.35 | 64.94 $\pm$ 0.28 | 71.77 $\pm$ 1.05 | 34.98 | 36.57
x |  | x |  |  66.72 $\pm$ 0.21 | 65.21 $\pm$ 0.28 | 65.00 $\pm$ 0.73 | 71.87 $\pm$ 1.14 | 35.22 | 36.50
x |  |  | x | 66.38 $\pm$ 0.23 | 65.83 $\pm$ 0.39 | 65.24 $\pm$ 0.70 | 70.96 $\pm$ 1.03 | 35.41 | 37.09

Comparing our results to the related work in Table 1 shows that our student
model with a performance of 37.74 clearly outperforms six out of eight methods
on ADE20K, even though it was trained using only the most basic distillation
loss. As the ADE20K dataset is the most complex dataset with 150 classes, this
observation is surprising.

##  6 Ablation Study

Another hyperparameter of KD is the initialization of student weights. The
effect of pre-training student backbone models on the ImageNet dataset
compared to random weight initialization has been studied with varying
conclusions. The authors of SKD find KD to be more efficient when the student
is initialized randomly [16]. CWD [26] contradicts this by stating that the
pre-trained weights help distillation, but calls the lift in student
performance less significant compared to the randomly initialized case. The
provided reason for this statement is that the relative improvement of the
student is smaller since the model was better when trained without a teacher.
[35] and [38] also report a higher absolute improvement for random
initialization. [30] is very unclear on the subject of weight initialization,
stating in a sketch of their algorithm that student weights are initialized
randomly, while the performance of student models trained without a teacher
suggests an initialization with pre-trained weights.

Our investigation of student weight initialization suggests that distillation
of teacher knowledge by the simple pixel-wise distillation loss does not
improve the performance of a randomly initialized student, while it does
improve the pre-trained one.

Table 7: mIoU of randomly initialized ResNet student trained on Cityscapes compared to literature. The best mIoU is highlighted in bold.  Method | SKD | CWD | Ours
---|---|---|---
Teacher | 78.56 | 78.5 | 78.24
$L_{CE}$ | 57.50 | 63.63 | 63.68
$L_{CE}+L_{PI}$ | 58.63 | - | 63.34
$L_{CE}+L_{PI}+L_{PA}+L_{HO}$ | 63.24 | 63.20 | -

Table 7 compares our findings to the results of [16] and [26]. Our experiments
find the validation mIoU of the student trained without a teacher to be 63.68,
averaged over four runs. This result is in line with CWD, but 6.18 percentage
points higher than what is reported by [16]. Also similar to CWD we find no
positive effect of distillation of teacher knowledge by means of the
investigated loss terms.

##  7 Environmental Impact

KD is an effective technique to reduce energy consumption at inference time,
as the student model requires less energy than the larger teacher. On the
other hand, energy consumption during student training is increased compared
to training the student without teacher. These two phenomena result in two
trade-offs: the first is performance loss compared to teacher vs. reduced
energy consumption during inference and the second is increased performance
compared to student-only vs. increased energy consumption during training. The
decision about when to use KD always depends on the exact use case. If a model
is expected to be deployed on a large number of devices, which all process
images at a high rate, energy use during training might be negligible compared
to inference and thus KD might be extremely beneficial.

Using the codecarbon.io tool we calculate the $CO_{2}$ emissions of training
and inference of our experiments with the ResNet student on Cityscapes and
show them in Table 8.

Table 8: $CO_{2}$ emissions of different model combinations at training and inference on Cityscapes. S is the ResNet student, numbers for inference are per 10,000 images. |  | training | inference
---|---|---|---
| mIoU | energy (kWh) |  $CO_{2}$ (kg) | energy (kWh) |  $CO_{2}$ (kg)
S | 70.45 | 0.841 | 0.298 | 0.161 | 0.057
S + KD | 72.75 | 2.52 | 0.891 | 0.161 | 0.057
T | 78.24 | 3.01 | 1.06 | 0.413 | 0.146

Training the student on its own takes 0.841kWh, which with the German
electricity conditions means 298g of emitted $CO_{2}$. As expected, the
addition of the teacher in KD leaves energy consumption during inference
unchanged, but increases energy consumption and emissions during training to
almost the three-fold. Comparing student and teacher reveals that the ResNet
student emits only 57g of $CO_{2}$ per 10,000 images at inference compared to
146g for the teacher model. This means, if we accept the decrease in
performance we could save over 60% of $CO_{2}$ emission during inference time.

##  8 Conclusion

In this work, we point out a significant comparability problem in the field of
KD for semantic segmentation, which will grow in relevance as more techniques
are published. We argue that comparibility can be improved by thoroughly
optimizing training hyperparameters and show that doing so eliminates the
gains from two accepted techniques, SKD and IFVD, on the two more complex of
the three investigated datasets. To facilitate an easier comparison in the
future we provide a detailed training protocol including optimal
hyperparameters for two student models and three datasets. As part of the
hyperparameter tuning we investigate the temperature parameter $\tau$, which
most publications in the field simply set to 1. We investigate the entropy of
class probability distributions of teacher output to visualize the softening
effect of $\tau>1$ and show that the temperature parameter can be beneficial
to the distillation of teacher knowledge.

Following up on this work and using the presented training with optimized
hyperparameters, it would be useful to provide a fair comparison of the other
loss terms in Table 1. Additionally, the understanding of KD in segmentation
could be deepened by an analysis of factors that decide when logit scaling is
helpful and why it improves distillation on Cityscapes and not on the other
datasets.

## Disclaimer

The results, opinions and conclusions expressed in this publication are not
necessarily those of Volkswagen Aktiengesellschaft.

## References

  * [1] An, S., Liao, Q., Lu, Z., Xue, J.H.: Efficient semantic segmentation via self-attention and self-distillation. IEEE Transactions on Intelligent Transportation Systems 23(9), 15256–15266 (2022). https://doi.org/10.1109/TITS.2021.3139001
  * [2] Badrinarayanan, V., Kendall, A., Cipolla, R.: Segnet: A deep convolutional encoder-decoder architecture for image segmentation. CoRR abs/1511.00561 (2015), http://arxiv.org/abs/1511.00561
  * [3] Beyer, L., Zhai, X., Royer, A., Markeeva, L., Anil, R., Kolesnikov, A.: Knowledge distillation: A good teacher is patient and consistent. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). pp. 10925–10934 (June 2022)
  * [4] Bucila, C., Caruana, R., Niculescu-Mizil, A.: Model compression. In: Proceedings of the 12th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. p. 535–541. KDD ’06, Association for Computing Machinery, New York, NY, USA (2006). https://doi.org/10.1145/1150402.1150464, https://doi.org/10.1145/1150402.1150464
  * [5] Cordts, M., Omran, M., Ramos, S., Rehfeld, T., Enzweiler, M., Benenson, R., Franke, U., Roth, S., Schiele, B.: The cityscapes dataset for semantic urban scene understanding. In: Proc. of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2016)
  * [6] Everingham, M., Gool, L.V., Williams, C.K.I., Winn, J.M., Zisserman, A.: The pascal visual object classes (voc) challenge. Int. J. Comput. Vis. 88(2), 303–338 (2010), http://dblp.uni-trier.de/db/journals/ijcv/ijcv88.html#EveringhamGWWZ10
  * [7] Feng, Y., Sun, X., Diao, W., Li, J., Gao, X.: Double similarity distillation for semantic image segmentation. IEEE Transactions on Image Processing 30, 5363–5376 (2021). https://doi.org/10.1109/tip.2021.3083113, https://doi.org/10.1109%2Ftip.2021.3083113
  * [8] Furlanello, T., Lipton, Z.C., Tschannen, M., Itti, L., Anandkumar, A.: Born again neural networks (2018). https://doi.org/10.48550/ARXIV.1805.04770, https://arxiv.org/abs/1805.04770
  * [9] Goodfellow, I.J., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y.: Generative adversarial networks (2014). https://doi.org/10.48550/ARXIV.1406.2661, https://arxiv.org/abs/1406.2661
  * [10] Hao, S., Zhou, Y., Guo, Y.: A brief survey on semantic segmentation with deep learning. Neurocomputing 406, 302–321 (2020). https://doi.org/https://doi.org/10.1016/j.neucom.2019.11.118, https://www.sciencedirect.com/science/article/pii/S0925231220305476
  * [11] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 770–778 (2016)
  * [12] He, T., Shen, C., Tian, Z., Gong, D., Sun, C., Yan, Y.: Knowledge adaptation for efficient semantic segmentation (2019). https://doi.org/10.48550/ARXIV.1903.04688, https://arxiv.org/abs/1903.04688
  * [13] Hinton, G., Vinyals, O., Dean, J.: Distilling the knowledge in a neural network (2015). https://doi.org/10.48550/ARXIV.1503.02531, https://arxiv.org/abs/1503.02531
  * [14] Ji, D., Wang, H., Tao, M., Huang, J., Hua, X., Lu, H.: Structural and statistical texture knowledge distillation for semantic segmentation. In: 2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). pp. 16855–16864 (2022). https://doi.org/10.1109/CVPR52688.2022.01637
  * [15] Li, H., Xiong, P., Fan, H., Sun, J.: Dfanet: Deep feature aggregation for real-time semantic segmentation. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) (June 2019)
  * [16] Liu, Y., Chen, K., Liu, C., Qin, Z., Luo, Z., Wang, J.: Structured knowledge distillation for semantic segmentation. In: 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). pp. 2599–2608 (2019). https://doi.org/10.1109/CVPR.2019.00271
  * [17] Long, J., Shelhamer, E., Darrell, T.: Fully convolutional networks for semantic segmentation. CoRR abs/1411.4038 (2014), http://arxiv.org/abs/1411.4038
  * [18] Menghani, G.: Efficient deep learning: A survey on making deep learning models smaller, faster, and better. CoRR abs/2106.08962 (2021), https://arxiv.org/abs/2106.08962
  * [19] Noh, H., Hong, S., Han, B.: Learning deconvolution network for semantic segmentation (2015). https://doi.org/10.48550/ARXIV.1505.04366, https://arxiv.org/abs/1505.04366
  * [20] Park, S., Heo, Y.S.: Knowledge distillation for semantic segmentation using channel and spatial correlations and adaptive cross entropy. Sensors 20(16) (2020). https://doi.org/10.3390/s20164616, https://www.mdpi.com/1424-8220/20/16/4616
  * [21] Paszke, A., Chaurasia, A., Kim, S., Culurciello, E.: Enet: A deep neural network architecture for real-time semantic segmentation. arXiv preprint arXiv:1606.02147 (2016)
  * [22] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., Chintala, S.: Pytorch: An imperative style, high-performance deep learning library. In: Advances in Neural Information Processing Systems 32, pp. 8024–8035. Curran Associates, Inc. (2019), http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf
  * [23] Romero, A., Ballas, N., Kahou, S.E., Chassang, A., Gatta, C., Bengio, Y.: Fitnets: Hints for thin deep nets (2015)
  * [24] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. CoRR abs/1505.04597 (2015), http://arxiv.org/abs/1505.04597
  * [25] Shannon, C.E.: A mathematical theory of communication. The Bell System Technical Journal 27(3), 379–423 (1948). https://doi.org/10.1002/j.1538-7305.1948.tb01338.x
  * [26] Shu, C., Liu, Y., Gao, J., Yan, Z., Shen, C.: Channel-wise knowledge distillation for dense prediction (2020). https://doi.org/10.48550/ARXIV.2011.13256, https://arxiv.org/abs/2011.13256
  * [27] Stanton, S., Izmailov, P., Kirichenko, P., Alemi, A.A., Wilson, A.G.: Does knowledge distillation really work? (2021). https://doi.org/10.48550/ARXIV.2106.05945, https://arxiv.org/abs/2106.05945
  * [28] Strudel, R., Garcia, R., Laptev, I., Schmid, C.: Segmenter: Transformer for semantic segmentation (2021). https://doi.org/10.48550/ARXIV.2105.05633, https://arxiv.org/abs/2105.05633
  * [29] Tan, M., Le, Q.: Efficientnet: Rethinking model scaling for convolutional neural networks. In: International conference on machine learning. pp. 6105–6114. PMLR (2019)
  * [30] Tian, Z., Chen, P., Lai, X., Jiang, L., Liu, S., Zhao, H., Yu, B., Yang, M.C., Jia, J.: Adaptive perspective distillation for semantic segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence 45(2), 1372–1387 (2023). https://doi.org/10.1109/TPAMI.2022.3159581
  * [31] Turc, I., Chang, M., Lee, K., Toutanova, K.: Well-read students learn better: The impact of student initialization on knowledge distillation. CoRR abs/1908.08962 (2019), http://arxiv.org/abs/1908.08962
  * [32] Wang, Y., Zhou, W., Jiang, T., Bai, X., Xu, Y.: Intra-class Feature Variation Distillation for Semantic Segmentation, pp. 346–362 (11 2020). https://doi.org/10.1007/978-3-030-58571-6_21
  * [33] Xie, E., Wang, W., Yu, Z., Anandkumar, A., Alvarez, J.M., Luo, P.: Segformer: Simple and efficient design for semantic segmentation with transformers (2021). https://doi.org/10.48550/ARXIV.2105.15203, https://arxiv.org/abs/2105.15203
  * [34] Xie, J., Shuai, B., Hu, J.F., Lin, J., Zheng, W.S.: Improving fast segmentation with teacher-student learning (2018). https://doi.org/10.48550/ARXIV.1810.08476, https://arxiv.org/abs/1810.08476
  * [35] Yang, C., Zhou, H., An, Z., Jiang, X., Xu, Y., Zhang, Q.: Cross-image relational knowledge distillation for semantic segmentation (2022). https://doi.org/10.48550/ARXIV.2204.06986, https://arxiv.org/abs/2204.06986
  * [36] Yang, Z., Li, Z., Shao, M., Shi, D., Yuan, Z., Yuan, C.: Masked generative distillation (2022). https://doi.org/10.48550/ARXIV.2205.01529, https://arxiv.org/abs/2205.01529
  * [37] Yu, C., Wang, J., Peng, C., Gao, C., Yu, G., Sang, N.: Bisenet: Bilateral segmentation network for real-time semantic segmentation. In: Proceedings of the European conference on computer vision (ECCV). pp. 325–341 (2018)
  * [38] Yuan, J., Qi, Q., Du, F., Intra, Z., Wang, F., Liu, Y.: Fakd: Feature augmented knowledge distillation for semantic segmentation (2022). https://doi.org/10.48550/ARXIV.2208.14143, https://arxiv.org/abs/2208.14143
  * [39] Yuan, L., Tay, F.E.H., Li, G., Wang, T., Feng, J.: Revisit knowledge distillation: a teacher-free framework. CoRR abs/1909.11723 (2019), http://arxiv.org/abs/1909.11723
  * [40] Zagoruyko, S., Komodakis, N.: Paying more attention to attention: Improving the performance of convolutional neural networks via attention transfer (2017)
  * [41] Zhang, Z., Zhou, C., Tu, Z.: Distilling inter-class distance for semantic segmentation (2022). https://doi.org/10.48550/ARXIV.2205.03650, https://arxiv.org/abs/2205.03650
  * [42] Zhao, H.: semseg. https://github.com/hszhao/semseg (2019), https://github.com/hszhao/semseg
  * [43] Zhao, H., Qi, X., Shen, X., Shi, J., Jia, J.: Icnet for real-time semantic segmentation on high-resolution images. In: Proceedings of the European conference on computer vision (ECCV). pp. 405–420 (2018)
  * [44] Zhao, H., Shi, J., Qi, X., Wang, X., Jia, J.: Pyramid scene parsing network. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 2881–2890 (2017)
  * [45] Zheng, S., Lu, J., Zhao, H., Zhu, X., Luo, Z., Wang, Y., Fu, Y., Feng, J., Xiang, T., Torr, P.H.S., Zhang, L.: Rethinking semantic segmentation from a sequence-to-sequence perspective with transformers (2020). https://doi.org/10.48550/ARXIV.2012.15840, https://arxiv.org/abs/2012.15840
  * [46] Zhou, B., Zhao, H., Puig, X., Fidler, S., Barriuso, A., Torralba, A.: Scene parsing through ade20k dataset. In: 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR). pp. 5122–5130 (2017). https://doi.org/10.1109/CVPR.2017.544

  *[mIoU]: Mean Intersection over Union


# Subspace Identification with Multiple Data Records: unlocking the archive

Chad M. Holcomb<EMAIL_ADDRESS>Robert R. Bitmead
<EMAIL_ADDRESS>Solar Turbines Inc., 4200 Ruffin Road, San Diego CA,
92123-5398, USA.  Department of Mechanical & Aerospace Engineering, University
of California, San Diego, 9500 Gilman Drive, La Jolla CA, 92093-0411, USA.

###### Abstract

We develop an approach to subspace system identification using multiple data
records and present a simple rank-based test for the adequacy of these data
for fitting the unique linear, noise-free, dynamic model of prescribed state-
vector, input-vector and output-vector dimensions. The approach is motivated
by the prospect of sorting through archives of operational data and extracting
a sequence of not-necessarily-contiguous data records individually
insufficient for providing identifiability but collectively making this
possible. The test of identifiability then becomes the sorting criterion for
accepting or rejecting new data records. En passant, the familiar Hankel
structure of the data matrices of subspace system identification is
reinterpreted and revised.

###### keywords:

Subspace; System Identification; Persistence of Excitation; Identifiabilty

††thanks: This paper was not presented at any IFAC meeting. Corresponding
author C. M. Holcomb.

, ,

##  1 Introduction

Dynamic System Identification is traditionally posed as the problem of fitting
the parameters of a dynamic model of a system using a single contiguous record
of input and output data. Because of the dynamics, breaking the record or
using two non-contiguous parts is disjunctive to the fitting process, since it
requires the treatment of additional initial conditions or the overwhelming of
initial condition effects via large data sets. We shall demonstrate that
segmented data need not stymie the identification. The vehicle for the
analysis is subspace system identification (SSI) in general (and the MOESP
algorithm in particular), which is amenable to multi-input/multi-output (MIMO)
system identification and which is increasingly the platform of choice for
many implementations of linear system identification.

The motivation for this problem stems from the advent of the cloud and the
industrial internet, where the penetration of remotely accessible data logging
capabilities and attendant expectations has expanded greatly. In contrast to
this new capability and the associated desire to use this copious data for
diagnostic and prognostic purposes lies the observation that in general
operational data is too poorly excited or too infrequently excited to reveal
the system dynamics reliably. Using a single very large data record swamps the
informative data with the under-informative. This feature has been examined by
Carrette et al. [1] with the conclusion that discarding uninformative data can
lead to improved estimation accuracy. By contrast, if archival data can be
used in pieces so as to yield collectively a well-excited experiment, then the
archive suddenly gains diagnostic interest and value. This is the driver for
the new methodology and for the associated test of adequacy of cumulative
excitation.

Peretzki et al. [8] provide a consistent analysis of the problems and values
of archived industrial data for the extraction of parametric system models.
They develop an approach to the assessment of signal data quality and
determination of its adequacy for model fitting by performing a running
condition number test on data segments. Their model structures consist of the
output regressing on Laguerre filtered input signals. Thus, the models do not
require initial states to accommodate non-contiguous segments. Their
motivation and data analytical focus parallel our approach here, although we
develop a subspace system identification algorithm and its attendant exact
theory for identifiability using ideas from behavioral systems theory in
comparison to their approach based on asymptotic normality.

Recently, other approaches to accommodate missing data in system
identification have been advanced. Markovsky [5, 4] considers small gaps
(including periodic gaps) in a single data record. The approach is to
subdivide the block-Hankel matrices of the single-record data equation into
complete block-Hankel sub-matrices and then to use linear algebraic approaches
to extract the underlying model. This yields a variant of SSI which focuses on
preserving the model fitting aspects. However, the excitation questions are
not broached.

The paper has two distinct parts. The first develops the precise theory by
dealing with multiple-segment system identification from noise-free data for
exact linear systems of known degree. Signal identifiability conditions are
presented in terms of rank. The second part presents an example from
industrial gas turbine systems, where the rank-revealing singular value
decomposition (SVD) replaces the rank calculations to accommodate noisy data
and approximate linearity.

Our approach is to present briefly the single-record SSI data equation and
approach in Section 2 before posing the multi-record formulation in Section 3,
drawing on their similarity. This multiple-segment algorithm is an extension
of that from [16, 15] associated with radio channel equalization. The
algorithm was presented in the form of this paper in [2] using the
instrumental variable variant, but without the current theory. Theorem 1 is
presented in Section 4 providing a testable sufficient identifiability
condition on the data in this noise-free, exact-modeling environment. These
results are a natural extension of those of Willems and Markovsky [11, 6] and
an improvement of the test of [2], which did not present a proof. A simple
computational example is given in this section to fix ideas. Then in Section
6, multiple archival gas turbine data sets are used to construct an
informative experiment from which a model can be calculated. These data sets
are not contiguous. Nor is any one data set completely informative by itself.

##  2 Single-record subspace system identification

As indicated above, we focus on the MOESP algorithmic approach to SSI [10] for
clarity. Since the underlying data equations, (9) and (31), are the same for
all SSI approaches, such as CVA [3], N4SID [7, 9], etc, the approach and data
adequacy test are applicable.

###  2.1 System

Consider a linear, time-invariant, noise-free, discrete-time, MIMO, dynamic
system $\mathcal{S}$, described by the state space system equations

$\displaystyle\mathcal{S}:=\left\\{\begin{array}[]{rl}x_{t+1}&=Ax_{t}+Bu_{t},\\\
y_{t}&=Cx_{t}+Du_{t},\end{array}\right.$ (3)

with input $u_{t}\in\mathbb{R}^{m}$, output $y_{t}\in\mathbb{R}^{p}$, and
state $x_{t}\in\mathbb{R}^{n}$.

######  Assumption 1.

For system $\mathcal{S}$ in (3) the dimensions $n$, $m$, $p$ are known and $n$
is the minimal state dimension. That is, $\mathcal{S}$ is observable and
controllable. Signals $\\{u_{t},y_{t}:t=1,\dots,N\\}$ are measured exactly.

Denote the following matrices.

$\displaystyle U_{\ell}$
$\displaystyle=\begin{bmatrix}u_{1}&u_{2}&u_{3}&\dots&u_{j}\\\
u_{2}&u_{3}&u_{4}&\iddots&u_{j+1}\\\ u_{3}&u_{4}&\iddots&&u_{j+2}\\\
\vdots&\iddots&&&\vdots\\\
u_{\ell}&u_{\ell+1}&u_{\ell+2}&\dots&u_{\ell+j-1}\end{bmatrix},$ (4)
$\displaystyle Y_{\ell}$
$\displaystyle=\begin{bmatrix}y_{1}&y_{2}&y_{3}&\dots&y_{j}\\\
y_{2}&y_{3}&y_{4}&\iddots&y_{j+1}\\\ y_{3}&y_{4}&\iddots&&y_{j+2}\\\
\vdots&\iddots&&&\vdots\\\
y_{\ell}&y_{\ell+1}&y_{\ell+2}&\dots&y_{\ell+j-1}\end{bmatrix},$ (5)

$\displaystyle\Gamma_{\ell}$ $\displaystyle=\begin{bmatrix}C\\\ CA\\\
\vdots\\\ CA^{\ell-1}\end{bmatrix},$ (6) $\displaystyle H_{\ell}$
$\displaystyle=\begin{bmatrix}D&0&0&0&\cdots&0\\\ CB&D&0&0&\cdots&0\\\
CAB&CB&D&0&\cdots&0\\\ \vdots&\vdots&\vdots&\ddots&&\vdots\\\
CA^{\ell-2}B&&&&\cdots&D\end{bmatrix},$ (7) $\displaystyle X_{\ell}$
$\displaystyle=\begin{bmatrix}x_{1}&x_{2}&\cdots&x_{j}\end{bmatrix}.$ (8)

Their dimensions are $\ell m\times j$, $\ell p\times j$, $\ell p\times n$,
$\ell p\times\ell m$, and $n\times j$, respectively. Note that: $U_{\ell}$ and
$Y_{\ell}$ are block-Hankel data matrices, $\Gamma_{\ell}$ is the extended
observability matrix, $H_{\ell}$ is a block-Toeplitz matrix of impulse
response parameters, and $X_{\ell}$ is a matrix of successive state values.

The noise-free data equation of subspace system identification using a single
segment is as follows.

$\displaystyle Y_{\ell}$
$\displaystyle=\Gamma_{\ell}X_{\ell}+H_{\ell}U_{\ell},$ (9)
$\displaystyle\begin{bmatrix}U_{\ell}\\\ Y_{\ell}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}I&0\\\
H_{\ell}&\Gamma_{\ell}\end{bmatrix}\begin{bmatrix}U_{\ell}\\\
X_{\ell}\end{bmatrix}.$ (10)

Subspace system identification proceeds in four steps from (9).

Step S.i: Orthogonal projection to $U_{\ell}.$ Define

$\displaystyle\Pi_{U_{\ell}}^{\perp}=I-U_{\ell}^{T}(U_{\ell}U_{\ell}^{T})^{{\dagger}}U_{\ell},$
(11)

and multiply the data equation (9) on the right to yield

$\displaystyle Y_{\ell}\Pi^{\perp}_{U_{\ell}}=$
$\displaystyle\Gamma_{\ell}X_{\ell}\Pi^{\perp}_{U_{\ell}}.$ (12)

The column space of $\Gamma_{\ell}$ is the column space of
$Y_{\ell}\Pi^{\perp}_{U_{\ell}}$.

Step S.ii: Extract the column space of $\Gamma_{\ell}$ using the singular
value decomposition.

$\displaystyle Y_{\ell}\Pi^{\perp}_{U_{\ell}}$
$\displaystyle=\begin{bmatrix}\mathcal{Q}_{1}&\mathcal{Q}_{2}\end{bmatrix}\begin{bmatrix}\Sigma_{1}&0\\\
0&\Sigma_{2}\end{bmatrix}\begin{bmatrix}\mathcal{P}_{1}^{T}\\\
\mathcal{P}_{2}^{T}\end{bmatrix},$ (13) $\displaystyle\hat{\Gamma}_{\ell}$
$\displaystyle=\mathcal{Q}_{1}\Sigma_{1}^{1/2}.$ (14)

[Note that weighting matrices can be included in this step which specifies
which variant of subspace system identification, N4SID, MOESP, CVA, etc, is
being used. This is immaterial to our analysis here.]

Step S.iii: Take the first $p$ rows of $\hat{\Gamma}_{\ell}$ to be $\hat{C}$.

$\displaystyle\hat{C}=\hat{\Gamma}_{l}(1:p,:),$ (15)

and solve a least-squares problem for $\hat{A}$.

$\displaystyle\hat{A}$
$\displaystyle=\left(\hat{\Gamma}_{\ell}(1:(\ell-1)p,:)\right)^{\dagger}\hat{\Gamma}(p+1:\ell
p,:).$ (16)

Step S.iv: Solve a least-squares problem for $\hat{B},$ $\hat{D},$ and
$\hat{x}_{1}$ as follows. Write

$\displaystyle y_{t}=$
$\displaystyle\sum\limits_{\tau=1}^{t-1}CA^{t-1-\tau}Bu_{\tau}+Du_{t}+CA^{t-1}x_{1}.$
(17)

Apply the vec operator and the Kronecker product $(\otimes)$ formula to (17)
to yield

$\displaystyle y_{t}$ $\displaystyle=\mbox{\rm vec}(y_{t}),$
$\displaystyle=\phi_{B}^{T}(t)\mbox{\rm vec}(B)+\phi_{D}^{T}(t)\mbox{\rm
vec}(D)+\phi_{x_{1}}^{T}(t)x_{1},$ (18)

where

$\displaystyle\phi_{B}^{T}(t)$
$\displaystyle=\sum_{\tau=1}^{t-1}u_{\tau}^{T}\otimes\hat{C}\hat{A}^{t-1-\tau},$
(19) $\displaystyle\phi_{D}^{T}(t)$ $\displaystyle=u_{t}^{T}\otimes I_{p},$
(20) $\displaystyle\phi_{x_{1}}^{T}(t)$ $\displaystyle=\hat{C}\hat{A}^{t-1}.$
(21)

Now assemble the output and regressor matrices as follows.

$\displaystyle\mathsf{Y}$
$\displaystyle=\begin{bmatrix}y_{1}^{T}&y_{2}^{T}&\dots&y_{N}^{T}\end{bmatrix}^{T},$
(22) $\displaystyle\Phi^{T}$
$\displaystyle=\begin{bmatrix}\Phi_{B}^{T}&\Phi_{D}^{T}&\Phi_{x_{1}}^{T}\end{bmatrix},$
(23)

with submatrices

$\displaystyle\Phi_{B}$
$\displaystyle=\begin{bmatrix}\phi_{B}(1)&\cdots&\phi_{B}(N)\end{bmatrix},$
(24) $\displaystyle\Phi_{D}$
$\displaystyle=\begin{bmatrix}\phi_{D}(1)&\cdots&\phi_{D}(N)\end{bmatrix},$
(25) $\displaystyle\Phi_{x_{1}}$
$\displaystyle=\begin{bmatrix}\phi_{x_{1}}(1)&\cdots&\phi_{x_{1}}(N)\end{bmatrix}.$
(26)

Define the parameter $(nm+mp+n)$-vector

$\displaystyle\theta$ $\displaystyle=\begin{bmatrix}\mbox{\rm
vec}(B)^{T}&\mbox{\rm vec}(D)^{T}&x_{1}^{T}\end{bmatrix}^{T}.$ (27)

The least-squares problem for $B$, $D$, and $x_{1}$ is

$\displaystyle\mathsf{Y}$ $\displaystyle=\Phi^{T}\theta,$

with estimate

$\displaystyle\hat{\theta}$ $\displaystyle=\begin{bmatrix}\mbox{\rm
vec}(\hat{B})\\\ \mbox{\rm vec}(\hat{D})\\\
\hat{x}_{1}\end{bmatrix}=(\Phi\Phi^{T})^{-1}\,\Phi\mathsf{Y}.$

##  3 Multi-record subspace system identification

In place of the single data record $\\{u_{t},y_{t}:1=1,2,\dots,N\\},$ of which
we use $\ell+j$ elements, suppose we have available $j$ data records of length
$\ell$:

$\bigg{\\{}\\{u_{t},y_{t}:t=t_{i}+1,2,\dots,t_{i}+\ell\\}:i=1,2,\dots,j\bigg{\\}}.$

Each record is a contiguous length-$\ell$ sequence of input and output signal
values from System $\mathcal{S}$. Each of the $j$ records has unknown initial
state vector $x_{t_{i}+1}$.

Define the following matrices

$\displaystyle\mathbb{U}_{\ell}$
$\displaystyle=\begin{bmatrix}u_{t_{1}+1}&u_{t_{2}+1}&u_{t_{3}+1}&\dots&u_{t_{j}+1}\\\
u_{t_{1}+2}&u_{t_{2}+2}&u_{t_{3}+2}&\dots&u_{t_{j}+2}\\\
u_{t_{1}+3}&u_{t_{2}+3}&u_{t_{3}+3}&\dots&u_{t_{j}+3}\\\
\vdots&\vdots&\vdots&&\vdots\\\
u_{t_{1}+\ell}&u_{t_{2}+\ell}&u_{t_{3}+\ell}&\dots&u_{t_{j}+\ell}\end{bmatrix},$
(28) $\displaystyle\mathbb{Y}_{\ell}$
$\displaystyle=\begin{bmatrix}y_{t_{1}+1}&y_{t_{2}+1}&y_{t_{3}+1}&\dots&y_{t_{j}+1}\\\
y_{t_{1}+2}&y_{t_{2}+2}&y_{t_{3}+2}&\dots&y_{t_{j}+2}\\\
y_{t_{1}+3}&y_{t_{2}+3}&y_{t_{3}+3}&\dots&{}_{t_{j}+3}\\\
\vdots&\vdots&\vdots&&\vdots\\\
y_{t_{1}+\ell}&y_{t_{2}+\ell}&y_{t_{3}+\ell}&\dots&y_{t_{j}+\ell}\end{bmatrix},$
(29) $\displaystyle\mathbb{X}_{\ell}$
$\displaystyle=\begin{bmatrix}x_{t_{1}+1}&x_{t_{2}+1}&\cdots&x_{t_{j}+1}\end{bmatrix}.$
(30)

Then a variant of (9) holds with the same $\Gamma_{\ell}$ and $H_{\ell}$. To
wit,

$\displaystyle\mathbb{Y}_{\ell}$
$\displaystyle=\Gamma_{\ell}\mathbb{X}_{\ell}+H_{\ell}\mathbb{U}_{\ell},$ (31)
$\displaystyle\begin{bmatrix}\mathbb{U}_{\ell}\\\
\mathbb{Y}_{\ell}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}I&0\\\
H_{\ell}&\Gamma_{\ell}\end{bmatrix}\begin{bmatrix}\mathbb{U}_{\ell}\\\
\mathbb{X}_{\ell}\end{bmatrix}.$ (32)

The central distinction between the single-record (9) and multi-record (31)
lies in the non-block-Hankel structure of the data matrices
$\mathbb{U}_{\ell}$ and $\mathbb{Y}_{\ell}$. Although, one could return to (9)
exactly from (31) by selecting data records with

$\displaystyle t_{1}=0,\;t_{2}=1,\;\dots,t_{j}=j-1.$

So the single-record analysis can be subsumed in the multi-record approach.

The solution algorithm for recovering $\hat{C}$ and $\hat{A}$ proceeds
precisely as in Steps S.(i-iii) of the single-record algorithm:

Step M.i: Multiply $\mathbb{Y}_{\ell}$ on the right by

$\displaystyle\Pi_{\mathbb{U}_{\ell}}^{\perp}=I-\mathbb{U}_{\ell}^{T}(\mathbb{U}_{\ell}\mathbb{U}_{\ell}^{T})^{{\dagger}}\mathbb{U}_{\ell}.$
(33)

Step M.ii Compute the SVD

$\displaystyle\mathbb{Y}_{\ell}\Pi_{\mathbb{U}_{\ell}}^{\perp}=\begin{bmatrix}\mathbb{Q}_{1}&\mathbb{Q}_{2}\end{bmatrix}\begin{bmatrix}\mathbb{S}_{1}&0\\\
0&\mathbb{S}_{2}\end{bmatrix}\begin{bmatrix}\mathbb{P}_{1}^{T}\\\
\mathbb{P}_{2}^{T}\end{bmatrix}$

and take

$\displaystyle\hat{\Gamma}_{\ell}$
$\displaystyle=\mathbb{Q}_{1}\mathbb{S}_{1}^{1/2}.$

Step M.iii: Extract $\hat{C}$ and $\hat{A}$ as earlier from
$\hat{\Gamma}_{\ell}$ via (15) and (16).

Step M.iv: Rewrite (17) for each record, $t=1,\dots,\ell,$ and record number,
$i=1,\dots,j,$

$\displaystyle y_{t_{i}+t}=$
$\displaystyle\sum\limits_{\tau=1}^{t-1}\hat{C}\hat{A}^{t-1-\tau}Bu_{t_{i}+\tau}+Du_{t_{i}+t}+\hat{C}\hat{A}^{t-1}x_{t_{i}+1}.$
(34)

Whence, using (19-21),

$\displaystyle y_{t_{i}+t}$ $\displaystyle=\phi_{B}^{T}(t_{i}+t)\mbox{\rm
vec}(B)+\phi_{D}^{T}(t_{i}+t)\mbox{\rm vec}(D)$ $\displaystyle\hskip
79.66771pt+\phi_{x_{t_{i}+1}}^{T}(t)x_{t_{i}+1},$ (35)

where we have deliberately indicated that the regressor terms for the initial
conditions depend only on $t$ and not on $t_{i}$.

Define an extended parameter $(nm+mp+nj)$-vector

$\displaystyle\vartheta$ $\displaystyle=\begin{bmatrix}\mbox{\rm
vec}(B)^{T}&\mbox{\rm
vec}(D)^{T}&x_{t_{1}+1}^{T}&x_{t_{2}+1}^{T}&\dots&x^{T}_{t_{j}+1}\end{bmatrix}^{T}.$
(36)

Then the least-squares problem for $B$, $D$ and the initial conditions has the
form

$\displaystyle\mbox{\rm vec}(\mathbb{Y}_{\ell})=\Upsilon\vartheta,$ (37)

with regression matrix, $\Upsilon$, of the structure

$\displaystyle\begin{bmatrix}\Phi_{B}^{T}(t_{1}+1)&\Phi_{D}^{T}(t_{1}+1)&\hat{C}&0&\dots&0\\\
\Phi_{B}^{T}(t_{1}+2)&\Phi_{D}^{T}(t_{1}+2)&\hat{C}\hat{A}&0&\dots&0\\\
\vdots&\vdots&\vdots&\vdots&&\vdots\\\
\Phi_{B}^{T}(t_{1}+\ell)&\Phi_{D}^{T}(t_{1}+\ell)&\hat{C}\hat{A}^{\ell-1}&0&\dots&0\\\
\Phi_{B}^{T}(t_{2}+1)&\Phi_{D}(t_{2}+1)&0&\hat{C}&\dots&0\\\
\vdots&\vdots&\vdots&\vdots&&\vdots\\\
\Phi_{B}^{T}(t_{2}+\ell)&\Phi_{D}^{T}(t_{2}+\ell)&0&\hat{C}\hat{A}^{\ell-1}&\dots&0\\\
\vdots&\vdots&\vdots&\vdots&&\vdots\\\
\Phi_{B}^{T}(t_{j}+1)&\Phi_{D}(t_{j}+1)&0&0&\dots&\hat{C}\\\
\vdots&\vdots&\vdots&\vdots&&\vdots\\\
\Phi_{B}^{T}(t_{j}+\ell)&\Phi_{D}^{T}(t_{j}+\ell)&0&0&\dots&\hat{C}\hat{A}^{\ell-1}\end{bmatrix}$
(38)

Then the least-squares estimate

$\displaystyle\hat{\vartheta}$
$\displaystyle=(\Upsilon^{T}\Upsilon)^{\dagger}\Upsilon^{T}\mbox{\rm
vec}(\mathbb{Y}_{\ell}),$ (39)

and $B$, $D$, $\\{x_{t_{i}+1}:i=1,\dots,j\\}$ are unpacked from
$\hat{\vartheta}$ using (36).

##  4 Identifiability from multiple data records

Our main result follows. Its proof is given in Section 5 as is the definition
of the maximal system lag, $L$, which is overbounded by the system state
dimension, $n$.

######  Theorem 1.

For $\ell>L,$ system $\mathcal{S}$’s maximal lag, the multiple record subspace
system identification algorithm applied to data satisfying Assumption 1 yields
a minimal state-space description of System $\mathcal{S}$ and all of its
corresponding initial state values, $\\{x_{t_{1}+1},\dots,x_{t_{j}+1}\\}$ if
the following two rank conditions are met.

$\displaystyle\text{\rm rank}\,\begin{bmatrix}\mathbb{U}_{\ell}\\\
\mathbb{Y}_{\ell}\end{bmatrix}$ $\displaystyle=m\ell+n,$ (40)

and

$\displaystyle\text{\rm rank}\,\mathbb{U}_{\ell}$ $\displaystyle=m\ell.$ (41)

We note that an immediate requirement of Theorem 1 is that the data matrices
$\mathbb{U}_{\ell}$ and $\mathbb{Y}_{\ell}$ contain sufficient columns,

$\displaystyle j\geq m\ell+n,$ (42)

that (40) and (41) might feasibly be satisfied. That is, at least $m\ell+n$
contiguous length-$\ell$ data records are required.

###  4.1 Noise-free computational example

To demonstrate the method’s capabilities with noise-free data and accurate
computation, we briefly present a computational example using MATLAB. We
select the discrete-time, single-input/single-output, noise-free, minimal,
second-order, linear time-invariant system with matrices

$\displaystyle A=\begin{bmatrix}0.9&0.2\\\ 0&0.8\end{bmatrix},\hskip
14.22636ptB=\begin{bmatrix}1\\\ 1\end{bmatrix},\hskip
14.22636ptC=\begin{bmatrix}1&1\end{bmatrix},\hskip 14.22636ptD=1.$

Here

$\displaystyle m=1,\hskip 14.22636ptn=2,\hskip 14.22636ptp=1.$

We select data matrix block dimension

$\displaystyle\ell=3>L=n=2.$

[We note that for a proper single-input/single-output system such as this, the
system maximal lag, $L,$ (developed in Section 5) equals the state dimension,
$n$, and thus $\ell=3$ is the minimal possible choice.]

Seven data sets are constructed from specified initial states, running for 20
time steps, with input signals each drawn from independent, uniform $[0,1]$
white noise processes. Data record selection vector

$\displaystyle\texttt{jayVec}=\begin{bmatrix}0&0&0&2&2&1&0\end{bmatrix},$

dictates that two length-3 data columns are drawn from Records 4 and 5, i.e.
the first four data from each, while one length-3 column is drawn from Record
6. The number of $(u_{t},y_{t})$ data pairs used is 11 and the number of
columns is

$\displaystyle j=5=m\ell+n.$

The initial states for these three records are chosen to be

$\displaystyle x_{4,1}=\begin{bmatrix}-1\\\ -1\end{bmatrix},\hskip
14.22636ptx_{5,1}=\begin{bmatrix}0.5\\\ 1\end{bmatrix},\hskip
14.22636ptx_{6,1}=\begin{bmatrix}1\\\ 0.5\end{bmatrix}.\hskip 14.22636pt$

The data matrix $\mathbb{U}_{\ell}$ is $3\times 5,$ has singular values
$(2.0142,0.2645,0.0906)$ and therefore is rank $m\ell=3$. Condition (41)
holds. Further, using svd to compute rank,

$\displaystyle\texttt{svd}\begin{bmatrix}\mathbb{U}_{\ell}\\\
\mathbb{Y}_{\ell}\end{bmatrix}=(10.6335,2.2668,0.2509,0.0998,0.0221),$

and condition (40), rank=5, holds for this data selection.

The identified state space realization matrices and initial states for each
record are

$\displaystyle\hat{C}$
$\displaystyle=\begin{bmatrix}-0.6084&0.7157\end{bmatrix},\hskip
14.22636pt\hat{A}=\begin{bmatrix}0.9428&1.0603\\\
-0.0058&0.7572\end{bmatrix},$ $\displaystyle\hat{B}$
$\displaystyle=\begin{bmatrix}-3.2885\\\ -0.0009\end{bmatrix},\hskip
14.22636pt\hat{D}=1.0000,$ $\displaystyle\hat{x}_{4,1}$
$\displaystyle=\begin{bmatrix}3.2885\\\
0.0009\end{bmatrix},\hat{x}_{5,1}=\begin{bmatrix}-2.5039\\\
-0.0326\end{bmatrix},\hat{x}_{6,1}=\begin{bmatrix}-2.4288\\\
0.0312\end{bmatrix}.$

We recover the state transformation matrix

$\displaystyle T=\begin{bmatrix}-0.3064&8.1890\\\
-0.3020&-7.4733\end{bmatrix},$

from which it is simply demonstrated that

$\displaystyle A=T\hat{A}T^{-1},\hskip 8.53581ptB=T\hat{B},\hskip
8.53581ptC=\hat{C}T^{-1},\hskip 8.53581pt\hat{D}=D,$
$\displaystyle\begin{bmatrix}x_{4,1}&x_{5,1}&x_{6,1}\end{bmatrix}=T\begin{bmatrix}\hat{x}_{4,1}&\hat{x}_{5,1}&\hat{x}_{6,1}\end{bmatrix}.$

Regression matrix, $\Upsilon,$ of (38) is $11\times 9$ with singular values
(3.0502, 1.7359, 1.3986, 1.1391, 1.0001, 1.0000, 0.6407, 0.2578, 0.0635). The
11 rows correspond to the data
$\\{y_{4,1},\,y_{4,2},\,y_{4,3},\,y_{4,4},\,y_{5,1},\,y_{5,2},\,y_{5,3},\,y_{5,4},\,y_{6,1},\,y_{6,2},\,y_{6,3}\\}$.
The 9 columns comprise 2 regressors for the elements of $\hat{B}$, one
regressor for $\hat{D}$, and 6 regressors for the elements of $\hat{x}_{4,1},$
$\hat{x}_{5,1},$ and $\hat{x}_{6,1}$.

##  5 Theory

We make the following definitions.

Matrix $W_{\ell,j}$:

For given $\ell$ and $j$,

$\displaystyle\hskip-14.22636ptw_{t}=\begin{bmatrix}u_{t}\\\
y_{t}\end{bmatrix}\;\;\text{and}\;\;\;W_{\ell,j}$
$\displaystyle=\begin{bmatrix}w_{1}&w_{2}&\dots&w_{j}\\\
w_{2}&w_{3}&\iddots&w_{j+1}\\\ \vdots&\iddots&&\vdots\\\
w_{\ell}&w_{\ell+1}&\dots&w_{\ell+j-1}\end{bmatrix}.$ (43)

Behavior of $\mathcal{S}$:

The set of all signal functions $\\{w_{t};t=1,\dots\\}$ satisfying
$\mathcal{S}$ is the behavior of $\mathcal{S}$, denoted $\mathfrak{B}$. The
length-$\ell$ behavior of system $\mathcal{S}$, denoted
$\mathfrak{B}_{[1,\ell]},$ is the span of the column space of $W_{\ell,j}$ as
initial states and signals $\\{w_{t}\\}$ vary over their entire range for all
signals satisfying (3).

Annihilator Polynomial Matrix $R(\sigma)$:

Minimal-degree, row-reduced $p\times(p+m)$ matrix polynomial $R(\sigma)$ in
backward shift operator $\sigma$ (unique up to left multiplication by a
unimodular polynomial matrix) satisfying

$\displaystyle R(\sigma)w_{t}=R^{\prime}(\sigma)\xi_{t},$ (44)

for all $w_{t}\in\mathfrak{B}$ and for auxiliary variables (states here)
$\xi_{t}$ and auxiliary $p\times n$ polynomial matrix $R^{\prime},$ is the
minimal annihilator polynomial of the behavior. Such $R(\sigma)$ yields the
left Matrix Fraction Description of the transfer function matrix of
$\mathcal{S}$ and, hence, uniquely identifies $\mathcal{S}$. This is known as
the kernel representation of $\mathcal{S}$. Thereom 1 of [12] establishes that
the presence of the states does not affect the choice of annihilator $R$.

Maximal lag $L$:

The maximal degree of the elements of $R(\sigma)$ is the maximal system lag,
$L$. Polynomial matrix $R(\sigma)$ has $L+1$ constant $p\times(p+m)$ matrix
coefficients. Evidently, $L\leq n,$ the system degree and minimal state
dimension.

We have the following property characterizing the dimension of the range of
$\mathfrak{B}_{[1,\ell]}$.

######  Lemma 1 (Willems (1986)).

For $\ell>L$ and minimal $\mathcal{S}$,

$\displaystyle{\rm dim}\,\mathfrak{B}_{[1,\ell]}=m\ell+n.$ (45)

Clearly the left nullspace of $W_{\ell,j}\supseteq{\rm
ker}\;\mathfrak{B}_{[1,\ell]}$ and, if rank $W_{\ell}={\rm
dim}\;\mathfrak{B}_{[1,\ell]},$ then left nullspace $W_{\ell,j}={\rm
ker}\;\mathfrak{B}_{[1,\ell]}$ and the kernel representation of $\mathcal{S}$
is recovered from this left nullspace. Whence,

######  Lemma 2 (Willems et al (2005)).

For $\ell>L$ and minimal $\mathcal{S}$ if,

$\displaystyle{\rm rank}\,W_{\ell,j}=m\ell+n,$ (46)

then the left nullspace of $W_{\ell,j}$ uniquely identifies system
$\mathcal{S}$.

The precise condition, $\ell>L$, which in [12, 6] is difficult to discern
because of somewhat mercurial notation, follows from the isomorphism between
the left nullspace of $W_{\ell,j}$ and the $\ell$ coefficients of $R(\sigma)$.
Since the degree of $R$ is $L$, there is requirement for at least $L+1$
coefficient matrices. So $\ell$ must equal or exceed $L+1.$

######  Lemma 3.

For $\ell>L$ and minimal $\mathcal{S}$, if (40),

$\displaystyle{\rm rank}\,\begin{bmatrix}\mathbb{U}_{\ell}\\\
\mathbb{Y}_{\ell}\end{bmatrix}$ $\displaystyle=m\ell+n,$

holds then the left nullspace of the multi-record data matrix
$\begin{bmatrix}\mathbb{U}_{\ell}\\\ \mathbb{Y}_{\ell}\end{bmatrix}$ uniquely
identifies $\mathcal{S}$.

Proof: The data matrix from (10) is a row permutation of the corresponding
matrix $W_{\ell,j}$. Thus the ranks and left nullities of the two matrices are
the same. Further, the columns of the multiple-record data matrix of (40) are
columns of $W_{\ell,j}$ for $j=t_{j}+\ell$ and this subset of columns spans
$\mathfrak{B}_{[1,\ell]}$ if condition (40) holds. Thus, the kernel
representation of $\mathcal{S}$ will then be recoverable from the left
nullspace of the multiple data matrix and the system is identifiable from this
data. $\square$

######  Lemma 4.

If minimal $\mathcal{S}$ is uniquely identifiable and $\ell\geq L$, then the
initial state for each column of $\begin{bmatrix}\mathbb{U}_{\ell}\\\
\mathbb{Y}_{\ell}\end{bmatrix}$ is uniquely calculable from that column.

Proof: This follows from the observability of $\mathcal{S}$ and the connection
between the maximal lag, $L$, and the observability indices [12, 14].
$\square$

### Proof of Theorem 1

From Lemma 3, (40) suffices for unique identifiability of $\mathcal{S}$ via
the nullspace of the multi-record data matrix. Condition (41) further ensures
that, within the data records, we are free to choose the $\\{u_{t}\\}$ as the
input signal. In the behavioral approach, Willems [12] defines the input
signals to be those elements which may be taken as free variables in the data;
a property associated with (41). More particularly, this condition assures
that $\Pi_{\mathbb{U}_{\ell}}^{\perp}$ of (33) has rank precisely $n$, thereby
permitting the exact recovery of the rank $n$ column space of the extended
observability matrix. $\square$

######  Remark 1.

We note that condition $\ell>L$ in Theorem 1 might be replaced by $\ell>n,$
the minimal state dimension, in case $L$ is unknown, which is frequently the
case in system identification. For SISO systems, $L=n$. But for MIMO systems,
$L\leq n$.

######  Remark 2.

It is instructive to compare the persistence of excitation condition derived
in [13] to (41). Willems et al. appeal to (46) in Lemma 2 for identifiability
and then seek a condition on the single-record
$\\{u_{t}:t=1,\dots,\ell+j-1\\}$ sequence alone which is capable of ensuring
(46). His condition is that

$\displaystyle{\rm rank}\;U_{\ell+n}=m(\ell+n),$ (47)

where $U_{\ell}$ is defined in (4). Interestingly, this latter condition
requires satisfaction of a rank property on signals, part of which do not
appear in the data set used for identification. One suspects that for quasi-
stationary data there is a frequency content interpretation.

From the perspective of multiple data records, a surprising alteration to the
theory is the disappearance of block-Hankel matrices. Single contiguous
length-$N$ data records yield $N-\ell+1$ candidate length-$\ell$ contiguous
sub-records which might be stacked into block-Hankel matrices – their column
order is unimportant to the problem. Yet Hankel matrices are fundamental to
much of Systems Theory and realizations. Property (47) preserves this focus on
Hankel.

######  Remark 3.

It is clear that the linear algebraic approach to the multi-record problem
parallels that of the single-record method with the inherent inclusion of more
initial conditions and a corresponding complication of the regression problem.
However, it is also evident from Lemma 4 that the estimation of these
additional initial states is no more problematic than the estimation of a
single initial state, since the corresponding regression matrix is a direct
sum of extended observability matrices, $\hat{\Gamma}_{\ell}$, which are full
column rank by construction. This mirrors the content of Thereom 1 of [12] in
showing that the kernel estimation problem is unencumbered by the presence of
initial states. In practice, the incorporation of additional initial
conditions can introduce numerical issues.

######  Remark 4 (The [low] price of multiple data segments).

For the noise-free exact computation of Section 4, the data segment selection
matrix `jayVec` was introduced to choose the number of columns of the data
matrix to be composed from each available segment – the example chooses two
columns from segments four and five plus a single column from segment six,
which involved precisely eleven input-output data pairs.

The example uses state dimension $n=2$ and minimal $\ell=3$ and minimal
$j=m\ell=5.$ Its second-order scalar transfer function with direct
feedthrough,

$\displaystyle P(z)$
$\displaystyle=\frac{b_{0}+b_{1}z^{-1}+b_{2}z^{-2}}{1+a_{1}z^{-1}+a_{2}z^{-2}},$

possesses exactly five independent parameters. By the same token, each initial
state requires two independent parameters.

Examining the exact minimal number of input-output data pairs implied by
(28-29) for this selection of $n$, $\ell$ and $j$, we see that: using a single
data segment requires exactly seven data; using two data segments requires
exactly nine data; three data segments requires eleven data, etc. So the
additional initial conditions implicitly force more data pairs to be used in
the construction of the data matrices, even though the data matrix dimension
does not alter.

##  6 Archival industrial gas turbine data sets

Part 2 of this paper applies these theoretically supported exact results based
on rank conditions to industrial data sets in order to demonstrate several
features:

  1. (i)

that the methods are amenable to handling less than pristine industrial data,
where linearity is only a local feature, models are admittedly approximate,
and data integrity is challenged;

  2. (ii)

that the rank conditions of the theory can be effectively replaced by rank-
revealing SVD conditions;

  3. (iii)

that the methods are applicable to MIMO systems directly; and

  4. (iv)

that model fitting can improve with the introduction of multiple data sets,
thereby unlocking the data archive.

The specific gas turbine system and its identified models are not examined in
detail since, at this stage, we are looking at the example as a proof of
concept for the methodology.

Data sets were collected from an operating 10 megawatt dry low-emissions
industrial gas turbine driving a centrifugal gas compressor in a natural gas
pipeline application. These data sets are historical and were not artificially
excited for system identification. In this case, data is sampled hourly and
aggregated remotely at a central monitoring facility spanning a two month
period. The gas turbine unit is periodically shut down for maintenance or due
to a system fault. Each shut-down event creates a gap in the data record. Of
the 17 contiguous data records, only 4 (records 5, 6, 15, 17) are of length
close to 200 hours: 245, 265, 172 and 193 respectively, thereby limiting both
the fidelity and model order achievable from a single data record. Each data
set is 2-input, fuel flow $W_{F}$ and bleed valve command $BV_{c}$, and
2-output, Stage 5 temperature $T_{5}$ and shaft speed $N_{GP}$. Individual
data sets 5 and 6 and then the multiple data sets $\\{5,6\\}$ and
$\\{5,6,15\\}$ were used to fit fourth-order models with

$\displaystyle\ell=5,\;\;m=2,\;\;n=4,\;\;p=2,$ (48)

using the algorithm of Section 3.

Figure 1 displays the singular values of the matrix $\mathbb{U}_{\ell}$ from
(28) while Figure 2 shows those of the combined input-output matrix
$\mathbb{W}_{\ell}=\begin{bmatrix}\mathbb{U}_{\ell}^{T}&\mathbb{Y}_{\ell}^{T}\end{bmatrix}^{T}$.
The minimal rank condition for $\mathbb{U}_{\ell}$ from (41) for the parameter
values (48) is $\ell m=10.$ The corresponding rank condition (40) for
$\mathbb{W}_{\ell}$ is $\ell m+n=14.$ These two figures demonstrate the
improvement in singular values achieved by incorporating further data sets
into the identification; the joint data sets possess uniformly larger singular
values than their constituents.

Figure 1: Singular value distributions of the input matrix $\mathbb{U}_{\ell}$
from (28) and Theorem 1, which requires the rank to be 10. Figure 2: Singular
value distributions of the input-output matrix $\mathbb{W}_{\ell}$ from
Theorem 1, which requires the minimal rank to be 14.

The dimensions of the regressor matrices $\Upsilon$ from (38) are:
490$\times$16 for the 245-data segment 5; 530$\times$16 for the 265-data
segment 6; and 1020$\times$20 for the combined data segments 5 and 6; and
1364$\times$24 for the combined $\\{5,6,15\\}$ segments. With the number of
rows given by number of input data (2 channels per time sample) and the number
of columns given by the parameters: 8 for $4\times 2$ $B$, 4 for $2\times 2$
$D$, and 4 for each $4\times 1$ initial state. In each case, $\Upsilon$ is
full rank computed by svd. Although for the multiple segments $\\{5,6,15\\}$
case, the identified $A$ has eigenvalues outside the unit circle, which leads
to ill-conditioning of $\Upsilon$.

The performance of the models identified using data set 5, data set 6, and
data sets $\\{5,6\\}$ is evaluated using data set 17, an independent
validation set of duration 193 hours. (The $\\{5,6,15\\}$ multi-segment model
was not considered because of its instability.) The MATLAB function compare
and data set 17 are used to evaluate the one-step-ahead prediction performance
of the three models. The resulting mean squared prediction error performance
for each output channel is presented in Table 1.

RMS error ($\times 10^{2}$) | Model5 | Model6 | Model{56}
---|---|---|---
$N_{GP}$ (output 1) | 1.1266 | 1.0895 | 1.0796
$T_{5}$ (output 2) | 1.0916 | 1.0980 | 0.9093
Table 1: Root mean square (RMS) one-step-ahead prediction error from the three
models with independent validation data set 17.

The conclusion to be reached from this analysis is that the capability to
accommodate multiple MIMO data segments into the system identification process
has a significant capacity to improve the quality of fit. In the specific gas
turbine example, this effect is evident in the improved data quality measured
by the singular values of the data matrices in Theorem 1 and largely due to
the increased number of input-output data pairs. By the same token, had one or
both data sets proven under-excited there is the possibility that the joint
experiment might be sufficiently excited.

##  7 Conclusion

A new formal approach to the incorporation of multiple data sets into linear
system identification has been developed and the behavioral theory of linear
systems applied and extended to yield easily tested sufficient conditions on
the collection of data sets for exact model identifiability. It is
demonstrated that the accommodation of the additional initial conditions
associated with non-contiguous data sets is achieved automatically by the
method. The motivation for this study is the availability of significant
quantities of industrial data in cloud-based archives of operational records,
which are likely fragmented and discontinuous in time while also being largely
under-excited. By advancing the prospect of system identification using
multiple data sets and by providing a sufficiency test for identifiability, we
are well placed to commence the profitable analysis of this archive. A proof-
of-concept MIMO example using industrial gas turbine data was provided and
demonstrated improved quantified model fitting performance.

Evidently, there remains considerable further work connected with the suite of
options associated with subspace system identification, such as the choice of
weighting matrices, instrumental variable steps, guaranteeing stability, etc.
Likewise, on the data management side, approaches are needed to questions of
automating the searching, sorting and selection of segments of data records to
achieve sufficiently informative sets from which to fit reliably models for
dognostics. The combination of industrial data archiving, multi-segment system
identification, and the computer science of data management should represent
an important foray into so-called Big Data.

## References

  * [1] P. Carrette, G. Bastin, Y.Y. Genin, and M. Gevers.  Discarding data may help in system identification.  IEEE Transactions on Signal Processing, 44:2300–2310, 1996.
  * [2] C.M. Holcomb and R.R. Bitmead.  Gas turbine analytics using archival data sets.  In Proce ASME Turbo Expo: GT2015, page to appear, Montreal, Canada, 2015.
  * [3] W.E. Larimore.  Canonical variate analysis in identification, filtering, and adaptive control.  In Proceedings of the 29th IEEE Conference on Decision and Control, volume 2, pages 596–604, 1990.
  * [4] I. Markovsky.  Approximate system identification with missing data.  In 52nd IEEE Conference on Decision and Control, pages 156–161, Florence Italy, 2013.
  * [5] I. Markovsky.  Exact system identification with missing data.  In 52nd IEEE Conference on Decision and Control, pages 151–155, Florence Italy, 2013.
  * [6] I. Markovsky, J.C. Willems, S. Van Huffel, and B. De Moor.  Exact and Approximate Modeling of Linear Systems: a Behavioral Approach.  Mathematical Modeling and Computation. SIAM, Philadelphia PA, 2006.
  * [7] P. Van Overschee and B. De Moor.  N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems.  Automatica, 30(1):75–93, 1994.
  * [8] D. Peretzki, A.J. Isaksson, A.C. Bittencourt, and K. Forsman.  Data mining of historic data for process identification.  In Proceedings of AIChE Annual Meeting, pages 1027–1033, Minneapolis MN, 2011.
  * [9] Peter Van Overschee and Bart De Moor.  Subspace identification for linear systems: theory, implementation, applications.  Kluwer Academic Publishers, Dordrecht, 1996.
  * [10] Michel Verhaegen.  Identification of the deterministic part of mmo state space models given in innovations form from input-output data.  Automatica, 30(1):61–74, Jan. 1994.
  * [11] J. Willems, P. Rapisarda, I. Markovsky, and B. De Moor.  A note on persistency of excitation.  Systems & Control Letters, 54(4):325–329, April 2005.
  * [12] Jan C. Willems.  From time series to linear system — Part I.. Finite dimensional linear time invariant systems.  Automatica, 22(5):561–580, 1986.
  * [13] J.C. Willems, P. Rapisarda, I. Markovsky, and B. De Moor.  A note on persistence of excitation.  Systems & Control Letters, 54(4):325–329, 2005.
  * [14] W.A. Wolovich.  Linear Multivariable Systems.  Springer-Verlag, New York, NY, 1974.
  * [15] Chengjin Zhang and R.R. Bitmead.  MIMO eqaulization with state-space channel models.  IEEE Transactions on Signal Processing, 56(10):5222–5231, 2008\.
  * [16] Chengjin Zhang and R.R.Bitmead.  Subspace system identification for training-based MIMO channel estimation.  Automatica, 41:1623–1632, 2005.