E-PANNs: Sound Recognition Using Efficient Pre-trained Audio Neural Networks

Arshdeep<EMAIL_ADDRESS>Haohe
<EMAIL_ADDRESS>Mark D. Plumbley333<EMAIL_ADDRESS>
Centre for Vision, Speech and Signal Processing (CVSSP),

University of Surrey, UK

ABSTRACT

Sounds carry an abundance of information about activities and events in our
everyday environment, such as traffic noise, road works, music, or people
talking. Recent machine learning methods, such as convolutional neural
networks (CNNs), have been shown to be able to automatically recognize sound
activities, a task known as audio tagging. One such method, pre-trained audio
neural networks (PANNs), provides a neural network which has been pre-trained
on over 500 sound classes from the publicly available AudioSet dataset, and
can be used as a baseline or starting point for other tasks. However, the
existing PANNs model has a high computational complexity and large storage
requirement. This could limit the potential for deploying PANNs on resource-
constrained devices, such as on-the-edge sound sensors, and could lead to high
energy consumption if many such devices were deployed. In this paper, we
reduce the computational complexity and memory requirement of the PANNs model
by taking a pruning approach to eliminate redundant parameters from the PANNs
model. The resulting Efficient PANNs (E-PANNs) model, which requires 36% less
computations and 70% less memory, also slightly improves the sound recognition
(audio tagging) performance. The code for the E-PANNs model has been released
under an open source license.

## . NTRODUCTION

Everyday sound environments include a wide range of sound activities and
events, such as traffic noise, road works, key jangling, music, coughing or
people talking. These environmental sound activities contain an abundance of
information and can potentially be used in various applications including
public security surveillance, monitoring activities in a home for assisted
living, healthcare, and improving the office, workplace and urban environment.

Recent advances in machine learning infrastructure and the availability of
large-scale dataset such as AudioSet [1] have attracted artificial
intelligence and machine learning (AI/ML) researchers to develop methods for
automatic sound activity recognition, commonly known as audio tagging. A
typical audio tagging system is shown in Figure 1. It takes audio recordings
from the surrounding using microphones and then, recognises various sound
activities that occur in the surrounding area. Audio tagging systems using
convolutional neural networks (CNNs) have shown promising performance compared
to the traditional hand-crafted methods [2]. However, CNNs are resource-hungry
due to the high computational cost arising from multiply-accumulate operations
(MACs) and from the memory requirement for CNNs. For example, one of the best
performing audio tagging networks from the pre-trained audio neural networks
(PANNs) [2] framework has approximately 81M parameters and requires more than
50G MACs for inference corresponding to a 10s audio clip. Due to this, it may
be challenging to deploy such large-scale CNNs on resource-constrained devices
having a limited power budget and limited memory, such as smart phones or
internet of things (IoT) devices. Moreover, when large-scale CNNs are used as
a feature extractor or as a classifier for other downstream tasks such as
acoustic scene classification [3], the high computational cost of the CNNs
makes them slow during inference and particularly during the training process,
where the large-scale CNNs may consume more energy and emit more
$\text{CO}_{2}$. For instance, an NVIDIA RTX-2080 Ti GPU used to train machine
learning models for 48 hours generates the equivalent amount of
$\text{CO}_{2}$ to that emitted by an average car driven for 13
miles444Machine learning $\text{CO}_{2}$ estimator:
https://mlco2.github.io/impact/##compute. Therefore, large-scale CNNs are not
efficient and environmental-friendly despite performing well.

Figure 1: An audio tagging system recognising various sound activities in a
given audio recording [4].

Our contributions: To reduce computational complexity and memory storage of
large-scale CNNs such as PANNs, we eliminate redundant parameters and reduce
computational complexity from the PANNs model by taking a pruning approach.
The resulting Efficient PANNs (E-PANNs) model has significantly less MACs and
memory storage with slightly improved performance compared to that of the
original PANNs model. We make a real-time sound recognition demonstration
using E-PANNs publicly available555https://github.com/Arshdeep-Singh-
Boparai/E-PANNs.git.

The rest of the paper is organised as follows. Section 2 introduces some
background including a brief overview of convolutional neural networks, and
background on existing audio tagging systems and the PANNs architecture.
Section 3 presents the proposed methodology used to obtain E-PANNs. Next, the
experimental setup and dataset used for experiments is explained in Section 4.
Section 5 presents experimental analysis. Finally, Section 6 concludes the
paper.

## . ackground

### 2.1. Convolutional Neural Networks (CNNs)

Convolutional neural networks (CNNs) are a type of artificial neural network
inspired by biological nervous systems such as the human brain. CNNs are
designed to learn from examples or from a dataset for an underlying task such
as classification of sound activity in the surrounding area. They learn
through an optimization process to update their parameters including weights
and filters across various types of intermediate layers, such as
convolutional, pooling, and dense layers. An architecture of a CNN is shown in
Figure 2. A convolutional layer has multiple feature maps, where each feature
map is produced by the convolutional operation on input and a filter.

In a CNN, filters are small matrices of size $(k\times k)$ with $c$ channels
that are convolved across the input data. The convolution operation, as given
in Equation 1, is a multiply-accumulate operation (MAC) that involves sliding
the filter, $\mathbf{F}$, over the input, performing element-wise
multiplications between the filter and the corresponding input patch $x$ of
size ($c\times k\times k$), and summing up the results to produce a single
value, $y$. This process is repeated for every possible position of the filter
across the input, resulting in a feature map. Similarly, other feature maps
are produced using other filters in a given convolutional layer. Subsequently
a bias $b$ and a non-linear activation function $f(.)$ is applied to the
elements of feature maps. Other than convolutional, pooling and dense layers,
some CNNs might have residual blocks, which contain skip or shortcut
connections that bypass one or more of the convolutional layers. For an
introduction to various architectures of CNNs, see the article [5].

Figure 2: Convolutional neural network (CNN) architecture comprising of
convolutional, pooling and dense layer.

$y=f(\sum_{i=1}^{c}\sum_{k_{1}=1}^{k}\sum_{k_{2}=1}^{k}(x_{i,k_{1},k_{2}}\times\mathbf{F}_{i,k_{1},k_{2}})+b).$
(1)

### 2.2. Existing Audio Tagging Frameworks

With the release of the large-scale AudioSet dataset [1], which has 2M audio
examples and 527 classes, several researchers have conducted studies to
improve the performance of neural networks on AudioSet classification and in
particular CNNs. Methods for AudioSet classification include CNNs [2], CNNs
with residual blocks [6] and and self-attention based Transformers [7]. A
summary of the existing methods alongwith their performance and the number of
parameters is given in Table 1. While Transformer-based methods perform better
than CNNs, Transformers have high computational complexity, which makes their
deployment on low-powered devices difficult compared to the CNNs. We shall
therefore focus on CNNs in this paper.

Table 1: Mean average precision (mAPs) obtained for AudioSet evaluation
dataset using various audio tagging systems with their number of parameters.

Method | Neural Network type | Parameters (x $10^{6}$) | mAPs
---|---|---|---
PANNs-CNN6 [2] | Plain | 4.8 | 0.343
PANNs-CNN10 [2] | Plain | 5.2 | 0.380
ERANNs-2-5 [6] | Residual | 38.2 | 0.446
ERANNS-1-6 [6] | Residual | 54.5 | 0.450
PANNs (ResNet38) [2] | Residual | 73.78 | 0.434
PANNs-CNN14 [2] | Plain | 80.75 | 0.431
PANNs (Wavegram-Logmel-CNN) [2] | Plain | 81.06 | 0.439
AST [8] | Transformer | 88.10 | 0.459
AST(ensemble) [8] | Transformer | 526.6 | 0.485

### 2.3. PANNs-CNN14 Architecture for AudioSet Classification

The primary motivation for the development of PANNs was to provide pre-
training systems for audio pattern recognition on extensive datasets, in this
case, the AudioSet dataset [1], which can be used as a baseline network to
extract features or for a classification task to other audio-related tasks.
The authors proposed several architectures for PANNs, including CNN-14, which
demonstrated promising performance on various audio pattern recognition tasks.
Experiments demonstrated that PANNs can generalize well to other tasks with
limited training data and outperform models trained from scratch on those
tasks.

The architecture of the PANNs-CNN14 model, including the number of parameters
across each convolutional layers, the total number of parameters and the model
size is shown in Figure 6 (a) in the Appendix. It has six convolutional
blocks, each with two convolutional layers, followed by a batch normalization
layer and a ReLU activation function. The number of filters in each
convolutional block gradually increases from 64 to 2048 from layer to layer.
Finally, there is a dense layer with a sigmoid activation function that
outputs the predicted probabilities for each class.

The PANNs-CNN14 model takes a log-mel spectrogram of size (1000 $\times$ 64)
computed from the 10-second-length audio input. The model is trained with data
augmentation techniques such as Mixup [9] and SpecAugment [10] for 600k
iterations. For more details on CNN14 and other PANNs models, see [2].

## . roposed Methodology to Reduce Complexity

To reduce the computational complexity and memory storage of CNNs such as
PANNs-CNN14, we use a filter pruning approach [11, 12, 13] that involves
removing or “pruning” unimportant filters from the network, i.e. these filters
that contribute less to the performance of the CNN. Filter pruning is inspired
by the idea that some filters in a CNN are redundant or have a negligible
contribution to the overall accuracy of the network [14, 15]. These filters
can be safely removed from the network without significantly affecting its
performance. For example, previous studies [16, 17] found that 73% of the
filters in a SoundNet network [18] do not provide discriminative information
across different acoustic scene classes, and eliminating such filters give
similar performance to that of using all the original filters in SoundNet. For
a survey on pruning techniques, see [11].

The typical steps in a filter pruning process are:

  1. 1.

For a baseline, train the original network to a desired level of accuracy, or
utilise an already existing pre-trained network;

  2. 2.

Rank the filters based on a certain “importance” criterion, such as their
relevance in contributing to the performance of the network;

  3. 3.

Remove the least important filters and their corresponding feature maps from
the network to obtain a pruned network;

  4. 4.

Fine-tune the pruned network to recover some of the lost accuracy.

There are several benefits of filter pruning in CNNs. Firstly, it can
significantly reduce the size of the network, reducing the memory footprint
and computation time required for inference by removing filters and the
corresponding feature maps generated by the filters. Secondly, the robustness
of the network in maintaining performance and generalization capabilities can
improve by removing filters that may be sensitive to small perturbations in
the input data.

In this paper, we will apply the filter pruning approach to the PANNs-CNN14
model, to produce an efficient PANNs model, E-PANNs. To obtain E-PANNs, we
follow a three step pipeline as shown in Figure 3. A detailed description of
the steps is given below,

(Step 1) Take a baseline PANNs network: We take the publicly available pre-
trained CNN14 from PANNs and denote it as PANNs-CNN14. PANNs-CNN14 has
approximately 81M parameters and requires 21G MACs666MACs computation Pytorch
package.: https://pypi.org/project/thop/ for inference of a 10s audio clip
with a performance metric, mAPs of 0.431.

(Step 2) Compute filter importance across convolutional layers of the baseline
network: We compute the relevance of the filters to decide whether to retain
or prune the filter for each layer of PANNs-CNN14 independently using (a),
(b), and (c) steps as shown in Figure 3.

Assume that there are $n$ filters in a given convolutional layer. With these
$n$ filters, first, we measure how well a given filter produces the output in
the convolutional layer. Our hypothesis is that a filter producing significant
output can capture specific patterns or features in the input significantly,
which would be useful for the subsequent convolutional layers of CNN to better
understand the input data at different levels of abstraction. Therefore, we
consider a filter producing significant output more important than others.
Then, we measure the importance of each filter and rank the filters according
to their importance in a given convolutional layer. Subsequently, the process
(a)-(c) is repeated for other convolutional layers as well. More information
about the importance calculation of the filters can be found in [19].

(Step 3) Obtaining E-PANNs: Once we obtain the importance of each filter in
various convolutional layers in PANNs-CNN14, we eliminate $p$ least important
filters and their corresponding feature maps from various convolutional
layers. This results in a pruned network. The pruned network is then re-
trained or fine-tuned with same data as used for training the baseline PANNs-
CNN14 network, to regain some of the lost performance owing to elimination of
the filters and their corresponding feature map. Finally, an efficient version
of PANNs-CNN14, E-PANNs is obtained.

Figure 3: A flow diagram to obtain E-PANNs.

## . xperimental Setup

In this section, we describe the experimental setup to prune the PANNs-CNN14
and then fine-tune the pruned network.

In PANNs-CNN14 model, the last six convolutional layers (C7 to C12) yield
approximately 99% of the total network parameters (See architecture in Figure
6 (a) in the Appendix), therefore we prune $p\in\\{25\%,50\%,75\%\\}$ filters
from the last six convolutional layers. The fine-tuning process of the pruned
PANNs-CNN14 is performed in the same way as that used to train the original
PANNs-CNN14 [2]. The optimization procedure uses the same loss function, data
augmentation and batch size etc, but with fewer iterations, which are
approximately 450k for the pruned PANNs-CNN14. To evaluate the performance of
the proposed E-PANNs, we test it on the audio tagging problem using the
AudioSet evaluation dataset.

## . erformance Analysis

Figure 4: Convergence during fine-tuning process of E-PANNs when different
number of filters are pruned from the baseline network. In brackets, the
maximum mAPs obtained is also given.

Figure 4 shows a convergence plot obtained during the fine-tuning process of
E-PANNs, when different numbers of filters are pruned from the baseline
network. We find that E-PANNs recovers the same performance as that of the
baseline network, even after pruning 25% and 50% filters, with only 100k and
200k iterations respectively. We also find that it uses at least 3 times fewer
iterations compared to the baseline network to achieve mAPs equal to 0.431.
Therefore, E-PANNs can be used as an alternative to the original PANNs with an
advantage of less computation and memory requirement.

Figure 5: mAPs and the number of parameters obtained using the proposed method
and the other existing CNN-based methods.

We find that pruning 25% of the filters across C7 to C12 layers of the
baseline network removes 41% of the parameters and 24% of the MACs with an
improved mAPs equal to 0.442, compared to that of the baseline network. Also
after pruning 50% of the filters, the mAPs obtained are better than that of
the baseline network, with 70% fewer parameters and 36% fewer MACs. The
architecture of the E-PANNs network obtained after pruning 50% of the filters
is shown in Figure 6 (b) in the Appendix. On the other hand, when 75% of the
filters are removed, the mAPs is 0.420 with 78% fewer parameters and 46% fewer
MACs. This suggests that there is a trade-off between the number of parameters
pruned and the mAPs. As shown in Figure 5 (see arrow trend), we find that when
the number of parameters in E-PANNs is less than 25% of the total number of
parameters in the baseline network, and that the mAPs obtained using E-PANNs
becomes less than the mAPs obtained using the baseline network. Therefore, the
maximum number of parameters that can be pruned from the baseline network,
PANNs-CNN14, without any degradation in performance are approximately 75%.

Figure 5 also shows the number of parameters and mAPs obtained using the other
existing CNN-based methods as given in Table 1. Our proposed method gives
better mAPs with 70% fewer parameters and 60% fewer parameters respectively
compared to the best performing PANNs-(Wavegram-Logmel-CNN) and the PANNs-
ResNet38 network. In general, compared to other methods such as PANNs-CNN10 or
ERANNs-1-6 [6], we find that there is a trade-off between the number of
parameters and the mAPs.

## . onclusion

Convolutional neural networks (CNNs) are effective for audio recognition, but
can be costly. We make these CNNs more efficient by reducing their
computational cost and memory storage. We present a framework to obtain
E-PANNs from a large-scale pre-trained audio neural network (PANNs) for audio
tagging. We find that removing few of the unimportant filters from original
PANNs reduces the computational complexity by 36% and 70% fewer parameters
with an improved performance. Therefore, E-PANNs can be used as an alternative
to the original PANNs for the downstream tasks for feature extraction or
classification due to the less computational complexity and low memory
requirement of E-PANNs compared to that of the PANNs.

## . cknowledgements

This work was partly supported by Engineering and Physical Sciences Research
Council (EPSRC) Grant EP/T019751/1 “AI for Sound (AI4S)”. For the purpose of
open access, the authors have applied a Creative Commons Attribution (CC BY)
licence to any Author Accepted Manuscript version arising.

## REFERENCES

  * 1. Jort F Gemmeke, Daniel PW Ellis, Dylan Freedman, Aren Jansen, Wade Lawrence, R Channing Moore, Manoj Plakal, and Marvin Ritter.  Audio set: An ontology and human-labeled dataset for audio events.  In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 776–780. IEEE, 2017.
  * 2. Qiuqiang Kong, Yin Cao, Turab Iqbal, Yuxuan Wang, Wenwu Wang, and Mark D Plumbley.  PANNs: Large-scale pretrained audio neural networks for audio pattern recognition.  IEEE/ACM Transactions on Audio, Speech, and Language Processing, 28:2880–2894, 2020.
  * 3. Irene Martín-Morató, Toni Heittola, Annamaria Mesaros, and Tuomas Virtanen.  Low-complexity acoustic scene classification for multi-device audio: Analysis of DCASE 2021 Challenge systems.  DCASE Workshop, 2021.
  * 4. Eduardo Fonseca, Manoj Plakal, Frederic Font, Daniel P. W. Ellis, and Xavier Serra.  Audio tagging with noisy labels and minimal supervision.  In DCASE2019 Workshop, NY, USA, 2019.
  * 5. Jiuxiang Gu, Zhenhua Wang, Jason Kuen, Lianyang Ma, Amir Shahroudy, Bing Shuai, Ting Liu, Xingxing Wang, Gang Wang, Jianfei Cai, et al.  Recent advances in convolutional neural networks.  In Pattern Recognition, volume 77, pages 354–377. Elsevier, 2018\.
  * 6. Sergey Verbitskiy, Vladimir Berikov, and Viacheslav Vyshegorodtsev.  ERANNS: Efficient residual audio neural networks for audio pattern recognition.  Pattern Recognition Letters, 161:38–44, 2022.
  * 7. Yuan Gong, Andrew Rouditchenko, Alexander H Liu, David Harwath, Leonid Karlinsky, Hilde Kuehne, and James R Glass.  Contrastive audio-visual masked autoencoder.  In International Conference on Learning Representations, 2023.
  * 8. Yuan Gong, Yu-An Chung, and James Glass.  AST: Audio Spectrogram Transformer.  In Proc. Interspeech 2021, pages 571–575, 2021.
  * 9. Hongyi Zhang, Moustapha Cisse, Yann N Dauphin, and David Lopez-Paz.  mixup: Beyond empirical risk minimization.  International Conference on Learning Representations (ICLR), 2017\.
  * 10. Daniel S Park, William Chan, Yu Zhang, Chung-Cheng Chiu, Barret Zoph, Ekin D Cubuk, and Quoc V Le.  SpecAugment: A simple data augmentation method for automatic speech recognition.  Interspeech, 2019.
  * 11. Tailin Liang, John Glossner, Lei Wang, Shaobo Shi, and Xiaotong Zhang.  Pruning and quantization for deep neural network acceleration: A survey.  Neurocomputing, 461:370–403, 2021.
  * 12. Jian-Hao Luo, Jianxin Wu, and Weiyao Lin.  ThiNet: A filter level pruning method for deep neural network compression.  In Proceedings of the IEEE International Conference on Computer Vision, pages 5058–5066, 2017.
  * 13. Haisong Ding, Kai Chen, and Qiang Huo.  Compressing CNN-DBLSTM models for OCR with teacher-student learning and tucker decomposition.  In Pattern Recognition, volume 96, page 106957. Elsevier, 2019.
  * 14. Jonathan Frankle and Michael Carbin.  The lottery ticket hypothesis: Finding sparse, trainable neural networks.  In International Conference on Learning Representations, 2019.
  * 15. Misha Denil, Babak Shakibi, Laurent Dinh, M.A. Ranzato, and Nando De Freitas.  Predicting parameters in deep learning.  In Advances in Neural Information Processing Systems, pages 2148–2156, 2013.
  * 16. Arshdeep Singh, Padmanabhan Rajan, and Arnav Bhavsar.  Deep hidden analysis: A statistical framework to prune feature maps.  In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 820–824. IEEE, 2019.
  * 17. Arshdeep Singh, Padmanabhan Rajan, and Arnav Bhavsar.  SVD-based redundancy removal in 1-D CNNs for acoustic scene classification.  In Pattern Recognition Letters, volume 131, pages 383–389. Elsevier, 2020.
  * 18. Yusuf Aytar, Carl Vondrick, and Antonio Torralba.  SoundNet: Learning sound representations from unlabeled video.  Advances in Neural Information Processing Systems, 29, 2016.
  * 19. Arshdeep Singh and Mark D Plumbley.  Efficient CNNs via Passive Filter Pruning.  arXiv preprint arXiv:2304.02319, 2023.

## ppendix

Architecture of baseline network: PANNs-CNN14 model and efficient PANNs-CNN14.

Figure 6: (a) The architecture of baseline network PANNs-CNN14. (b) The
architecture of E-PANNs obtained after removing 50% of the filters from C7 to
C12 convolutional layers of PANNs-CNN14. The number of parameters for each
convolutional layer and the parameter (model) size is also shown. The PANNs-
CNN14 has 12 trainable convolutional layers, denoted as C1 to C12. The first
two convolutional layers (conv1d-1. conv1d-2) are non-trainable and their
purpose is to convert the raw audio signal to log-mel spectrogram.


# A Dynamic Observer for a Class of Infinite-Dimensional Vibrating Flexible
Structures

Alexander Zuyev1,2,3 and Julia Kalosha2 1Max Planck Institute for Dynamics of
Complex Technical Systems, Sandtorstraße 1, 39106 Magdeburg, Germany
<EMAIL_ADDRESS>of Applied Mathematics & Mechanics,
National Academy of Sciences of Ukraine, G. Batiuka St., 19, Sloviansk,
Ukraine 84100<EMAIL_ADDRESS>von Guericke University Magdeburg,
Universitätsplatz 2, 39106 Magdeburg, Germany
This work was supported in part by by the German Research Foundation (DFG)
under Grant ZU 359/2-1.

###### Abstract

Infinite-dimensional control systems with outputs are considered in the
Hamiltonian formulation with generalized coordinates. An explicit scheme for
constructing a dynamic observer for this class of systems is proposed with
arbitrary gain coefficients. Sufficient conditions for the convergence of the
constructed observer are obtained on the basis of the invariance principle.
This result is applied to a flexible beam model attached to a mass-spring
system with lumped and distributed actuators. The estimation error decay is
illustrated with numerical simulations of finite-dimensional approximations of
the observer dynamics.

##  I INTRODUCTION

Since the full state of a dynamical system is not available for direct
measurement in many applied problems, it becomes crucial for optimal control
and stabilization to estimate the state in a reliable manner employing
observations of a limited number of system outputs. A classical method of
estimating the state vector using the inputs and outputs, proposed by D.
Luenberger [1], has been widely extended for various classes of control
systems, including distributed parameter systems [2]–[4]. Note that the Kalman
observability rank condition or Hautus test are not directly applicable for
control systems with infinite degrees of freedom. So, powerful operator
techniques have been developed to unfold infinite-dimensional modifications of
classical methods (see, e.g.,[5]– [9]).

In recent years, a variety of results in observer design and observer-based
control have been obtained for specific engineering problems. A Luenberger-
type observe is proposed in [10] with the use of the Lyapunov method for a
flexible structure modeled by the Timoshenko beam. A feedback control scheme
based on the disturbance observer is designed for a flexible spacecraft
subject to external disturbance, input magnitude, and rate constraints in
[11]. A reduced-order nonlinear state observer is proposed in [12] for a
flexible-link multibody system by using model reduction techniques. In the
paper [13], a nonlinear observer based control is presented for a rotating
flexible structure modeled as the Timoshenko beam. The exponential convergence
of the estimation scheme and the asymptotic stability of the closed-loop
system are proved in this work. A state estimator for the motion of a slewing
flexible structure with rigid body is designed in [14]. In [15], the
Luenberger observer is used for estimating the internal state of a “black-box”
linear system, and the optimal feedback control is proposed. Stability under
the Luenberger observer parametrization is proved based Lyapunov stability
theory. In all these works results are experimentally validated.

The observability and controllability of diagonal systems with the infinite-
dimensional state and finite-dimensional output are investigated in [16],
[17]. In these papers, necessary and sufficient conditions for the exact
observability are extended to the systems that are not necessarily
exponentially stable. The above results are related to an infinite-dimensional
version of the Hautus test proposed by D.L. Russell and G. Weiss in [18]. Note
that in the works of B. Jacob and H. Zwart, the exact observability problem is
solved for diagonal systems without input and with the output signal of class
$L^{2}(0,\infty)$.

In the present work, we will propose an explicit observer design for a
controlled Hamiltonian system without the requirement that the output signal
is in $L^{2}(0,\infty)$. While our assumptions are weaker than in [16], [17],
we do not formulate exact observability conditions and prove that the
observation error tends to zero only asymptotically.

This paper may be contemplated as a further extension of the approach
developed in [19] to distributed parameter systems. In Section II, we present
a linear control system with the parameterized Hamiltonian operator and a
finite-dimensional output. It is shown that the considered system generates a
$C_{0}$-semigroup of operators in the real Hilbert space ${\ell^{2}}$ under a
suitable assumption on the parameters. Then we describe the observer design
procedure in Section III. Our stability analysis of the error dynamics is
based on LaSalle’s invariance principle (or an infinite-dimensional version of
the Barbashin–Krasovskii theorem) in Section V. A key part in the stability
proof (proof of Theorem 1) requires the precompactness of trajectories, which
is formally addressed in Section IV. Conditions for the observer convergence
are discussed in terms of parameters of the control system.

Section VI is devoted to a particular case, if the considered class of systems
describes vibrations of the flexible beam with an attached mass discussed in
[20]. In this case, parameters of the infinitesimal generator arise as
eigenvalues of a fourth-order differential operator. It has been proved in
[19] that the above operator is self-adjoint and positive, so its eigenvalues
form a sequence of positive real numbers. Besides, it has been proved in [21]
that all roots of the corresponding characteristic equation are simple and
satisfy certain growth condition. As a consequence of the eigenvalue
distribution added to a natural assumption that at least one output
coefficient is nonzero for each mode, the observer convergence is established
in Theorem 2. The observation error dynamics for the considered beam model is
illustrated by numerical simulations in Section VII.

##  II INFINITE-DIMENSIONAL CONTROL SYSTEM WITH OUTPUT

Consider an infinite-dimensional control system of the form

$\dot{z}=Az+Bu,\quad z\in H=\ell^{2}\times\ell^{2},\quad
u\in{\mathbb{R}}^{k+1},$ (1) $y=Cz,\quad y\in{\mathbb{R}}^{r},\quad
k,r\in\mathbb{N},$ (2)

where the state vector $z=\left(\begin{array}[]{c}\xi\\\ \eta\\\
\end{array}\right)$ consists of two components
${\xi=(\xi_{1},\xi_{2},\dots)^{T}\in\ell^{2}}$ and
${\eta=(\eta_{1},\eta_{2},\dots)^{T}\in\ell^{2}}$, the control
$u=(u_{0},u_{1},\dots,u_{k})^{T}$ and the output $y=(y_{1},\dots,y_{r})^{T}$
are finite-dimensional vectors. The inner product in the real Hilbert space
$H$ is inherited from $\ell^{2}$:

$\left\langle\left(\begin{array}[]{c}\xi\\\ \eta\\\
\end{array}\right),\left(\begin{array}[]{c}\tilde{\xi}\\\ \tilde{\eta}\\\
\end{array}\right)\right\rangle_{H}=\sum_{j=1}^{\infty}\left(\xi_{j}\tilde{\xi}_{j}+\eta_{j}\tilde{\eta}_{j}\right),$

and the operators $A:D(A)\to H$, $B:\mathbb{R}^{k+1}\to H$,
$C:H\to\mathbb{R}^{r}$ are defined in terms of the $\xi$\- and
$\eta$-component of $z\in H$ as

$A:z=\begin{pmatrix}\xi\\\ \eta\end{pmatrix}\mapsto
Az=\begin{pmatrix}\Omega\eta\\\
-\Omega\xi\end{pmatrix},\;Bu=\begin{pmatrix}0\\\
B_{1}u\end{pmatrix},\;Cz=C_{1}\xi,$

where $\Omega={\rm diag}(\omega_{1},\omega_{2},...),$ and the operators
$B_{1}:\mathbb{R}^{k+1}\to{\ell^{2}}$, $C_{1}{\ell^{2}}\to\mathbb{R}^{r}$ are
given by their matrices

$B_{1}=\left(\begin{array}[]{cccc}b_{10}&b_{11}&\ldots&b_{1k}\\\
b_{20}&b_{21}&\ldots&b_{2k}\\\
\vdots&\vdots&\vdots&\vdots\end{array}\right),\quad\sum\limits_{j=1}^{\infty}b_{ji}^{2}<\infty,\;i=0,\dots,k,$
$C_{1}=\left(\begin{array}[]{llc}c_{11}&c_{12}&\ldots\\\
\vdots&\vdots&\vdots\\\
c_{r1}&c_{r2}&\ldots\end{array}\right),\quad\sum\limits_{j=1}^{\infty}c_{sj}^{2}<\infty,\;s=1,\dots,r.$

Thus,

$D(A)=\left\\{z=\left(\begin{array}[]{c}\xi\\\ \eta\\\ \end{array}\right)\in
H:\sum\limits_{j=1}^{\infty}\omega_{j}^{2}(\xi_{j}^{2}+\eta_{j}^{2})<\infty\right\\}.$

System (1) is considered as a mathematical model of vibrating flexible
structures with force actuation and displacement measurement. We assume that
the spectrum of $A$ is purely imaginary without resonances, which is
formalized in the following way.

###### Assumption 1

The diagonal entries $\omega_{j}$ of $\Omega$ are mutually distinct positive
real numbers such that

${0<\omega_{1}<\omega_{2}<\dots<\omega_{n}<\dots}\;.$

It will be shown in section VI that Assumption 1 is fulfilled for the
particular system due to properties of the infinitesimal generator.

Unlike finite-dimensional ODE systems with smooth right-hand side, for which
the conditions of existence and uniqueness of the Cauchy problem solution are
clear, the question about well-posedness of the Cauchy problem for infinite-
dimensional systems becomes one of the major tasks. For dealing with this the
semigroup representation is commonly used.

###### Lemma 1

The operator $A$ generates a $C_{0}$-semigroup $\\{e^{tA}\\}_{t\geq 0}$ on
$H$.

###### Proof:

The domain $D(A)$ is dense in $H$. For any $z\in D(A)$, we have

$\langle Az,z\rangle_{H}\equiv 0,$ (3)

so the map $A$ is dissipative.

The inverse operator $A^{-1}:H\to H$ is defined as follows:

$A^{-1}=\left(\begin{array}[]{cc}0&-\tilde{A}_{1}\\\ \tilde{A}_{1}&0\\\
\end{array}\right),$

where $\tilde{A}_{1}={\rm
diag}\left(\frac{1}{\omega_{1}},\frac{1}{\omega_{2}},\dots\right)$.

It is easy to verify that the following estimate is fulfilled:

$\|A^{-1}z\|^{2}\leq\frac{1}{\omega_{1}^{2}}\;\|z\|^{2},\quad\forall\;z\in H,$

so $A^{-1}$ is bounded, therefore the operator $A$ is closed in $H$.

The resolvent $Res(A)$ of $A$ may be constructed as the map $Res(A)=(I-\lambda
A)^{-1}:H\to H$ with the components

$Res(A)=\left(\begin{array}[]{cc}R_{1}&R_{2}\\\ -R_{2}&R_{1}\\\
\end{array}\right),$

where

$\displaystyle R_{1}$ $\displaystyle={\rm
diag}\left(\frac{1}{\lambda^{2}\omega_{j}^{2}+1},\;j=1,2,\dots\right),$
$\displaystyle R_{2}$ $\displaystyle={\rm
diag}\left(\frac{\lambda\omega_{j}}{\lambda^{2}\omega_{j}^{2}+1},\;j=1,2,\dots\right)$

for some $\lambda>0$. Here and in the sequel, $I$ is the identity operator on
$H$.

As the transformation $z\mapsto Res(A)(z)$ is defined for any $z\in H$, the
range of $I-\lambda A$ coincides with $H$. Thus, the map $A$ is maximal.
Moreover, because of condition (3), $A$ is $m$-dissipative. Being densely
defined, $m$-dissipative, and closed, the operator $A$ satisfies the
conditions of the Lumer–Phillips theorem and thus generates a
$C_{0}$-semigroup on $H$. ∎

##  III LUENBERGER OBSERVER DESIGN

Our aim is to construct a Luenberger-type observer for system (1) in the form

$\dot{\bar{z}}(t)=(A-FC)\bar{z}(t)+Bu(t)+Fy(t),$ (4)

such that, for any initial conditions $z(0),\bar{z}(0)\in H$ and any
admissible control $u:[0,+\infty)\to{\mathbb{R}}^{k+1}$, the corresponding
solutions $z(t)$ and $\bar{z}(t)$ of (1), (2) and (4) satisfy the property

$\|z(t)-\bar{z}(t)\|\to 0\quad\;\text{as}\;t\to+\infty.$ (5)

For this purpose, we will define the operator $F:{\mathbb{R}}^{r}\to H$ as
$F=\left(\begin{array}[]{c}f\\\ g\\\ \end{array}\right)$ with
$f,g:{\mathbb{R}}^{r}\to{\ell^{2}}$ given by their matrices

$f=\left(\begin{array}[]{ccc}f_{11}&\dots&f_{1r}\\\ f_{21}&\dots&f_{2r}\\\
\vdots&\vdots&\vdots\\\ \end{array}\right)\quad\text{and}\quad
g=\left(\begin{array}[]{ccc}g_{11}&\dots&g_{1r}\\\ g_{21}&\dots&g_{2r}\\\
\vdots&\vdots&\vdots\\\ \end{array}\right),$

where

$f_{js}=\gamma_{s}c_{sj},\quad g_{js}=0,\quad s=\overline{1,r},\quad
j=1,2,\dots,$ (6)

and $\gamma_{s}>0$ are gain parameters.

Let us consider the observation error ${e(t)=z(t)-\bar{z}(t)}$ and denote
$e(t)=\left(\begin{array}[]{c}\Delta\\\ \delta\\\ \end{array}\right)$, where
${\Delta=(\Delta_{1},\Delta_{2},\dots)^{T}}$,
${\delta=(\delta_{1},\delta_{2},\dots)^{T}}$. Then the error dynamics can be
written as follows:

$\dot{e}(t)=(A-FC)e(t).$ (7)

The property (5) is equivalent to

$\|e(t)\|\to 0\quad\text{as}\quad t\to\infty,$

where $e(t)$ is the solution of (7).

##  IV PRECOMPACTNESS OF THE TRAJECTORIES

In this section, we will construct the resolvent $Res(\hat{A})$ of the
operator $\hat{A}=A-FC$ and prove its compactness for some $\lambda>0$. This
result will be then applied in Section V for the proof of our main result
(Theorem 1).

In order to construct $Res(\hat{A})$, we solve the equation

$(\hat{A}-\lambda I)\left(\begin{array}[]{c}\Delta\\\ \delta\\\
\end{array}\right)=\left(\begin{array}[]{c}\bar{\Delta}\\\ \bar{\delta}\\\
\end{array}\right)$ (8)

with respect to $(\Delta,\delta)^{T}\in H$ for any given
$(\bar{\Delta},\bar{\delta})^{T}\in H$, $\lambda>0$.

Let us define the matrix

$M=(M_{sp})=\lambda\gamma_{p}\sum\limits_{i=1}^{\infty}\frac{c_{si}c_{pi}}{\lambda^{2}+\omega_{i}^{2}}+{\bf\delta_{sp}},\quad
s,p=\overline{1,r}$

and its inverse $M^{-1}=(M^{-1}_{sp})_{s,p=\overline{1,r}}$. The symbol
${\bf\delta_{sp}}$ denotes the Kronecker delta.

Note that, for small $\lambda>0$, the matrix $M$ is a small disturbance of the
identity matrix, so the inverse matrix $M^{-1}$ exists and satisfies the norm
estimate

$\|M^{-1}\|\leq 1+\mathcal{O}(\lambda)\quad\text{as}\;\lambda\to 0.$

The solution of (8) can be presented as follows:

$\displaystyle\Delta_{j}$
$\displaystyle=-\frac{1}{\lambda^{2}+\omega_{j}^{2}}\left(\lambda\bar{\Delta}_{j}+\omega_{j}\bar{\delta}_{j}+\lambda\sum\limits_{s=1}^{r}\gamma_{s}c_{sj}\phi_{s}\right),$
(9) $\displaystyle\delta_{j}$
$\displaystyle=\frac{1}{\lambda^{2}+\omega_{j}^{2}}\left(\omega_{j}\bar{\Delta}_{j}-\lambda\bar{\delta}_{j}+\omega_{j}\sum\limits_{s=1}^{r}\gamma_{s}c_{sj}\phi_{s}\right).$

The parameters $\phi_{s}$, $s=\overline{1,r}$, can be obtained from the
following equation:

$\left(\phi_{1},\dots,\phi_{r}\right)^{T}=M^{-1}\bar{M},$

where $\bar{M}=-{\rm
colon}\left(\sum\limits_{i=1}^{\infty}\frac{c_{si}}{\lambda^{2}+\omega_{i}^{2}}\left(\lambda\bar{\Delta}_{i}+\omega_{i}\bar{\delta}_{i}\right),\;s=\overline{1,r}\right)$.
So,

$\phi_{s}=-\sum\limits_{p=1}^{r}M^{-1}_{sp}\sum\limits_{i=1}^{\infty}\frac{c_{pi}}{\lambda^{2}+\omega_{i}^{2}}\left(\lambda\bar{\Delta}_{i}+\omega_{i}\bar{\delta}_{i}\right)\quad
s=\overline{1,r}.$

The resolvent of $\hat{A}$ is acting as

$Res(\hat{A})=\left(\begin{array}[]{cc}R^{1}&R^{2}\\\ R^{3}&R^{4}\\\
\end{array}\right),$

where

$\displaystyle
R^{1}_{ji}=\frac{\lambda}{\lambda^{2}+\omega_{j}^{2}}\left(\frac{\lambda}{\lambda^{2}+\omega_{i}^{2}}\sum\limits_{s,p=1}^{r}\gamma_{s}M^{-1}_{sp}c_{sj}c_{pi}-{\bf\delta_{ji}}\right),$
$\displaystyle
R^{2}_{ji}=\frac{1}{\lambda^{2}+\omega_{j}^{2}}\left(\frac{\lambda\omega_{i}}{\lambda^{2}+\omega_{j}^{2}}\sum\limits_{s,p=1}^{r}\gamma_{s}M^{-1}_{sp}c_{sj}c_{pi}-\omega_{j}{\bf\delta_{ji}}\right),$
$\displaystyle R^{3}_{ji}=-\frac{\omega_{j}}{\lambda}R^{1}_{ji},$
$\displaystyle
R^{4}_{ji}=-\frac{1}{\lambda^{2}+\omega_{j}^{2}}\left(\frac{\omega_{j}\omega_{i}}{\lambda^{2}+\omega_{i}^{2}}\sum\limits_{s,p=1}^{r}\gamma_{s}M^{-1}_{sp}c_{sj}c_{pi}+\lambda{\bf\delta_{ji}}\right).$

As Assumption 1 is not be sufficient to guarantee the convergence of series
appearing in (9), we introduce the following

###### Assumption 2

The series $\sum\limits_{i=1}^{\infty}\frac{1}{\omega_{i}^{2}}$ is convergent.

###### Lemma 2

If Assumption 2 is fulfilled, then each positive semitrajectory
$\\{e(t)\\}_{t\geq 0}$ of system (7) is precompact in $H$.

###### Proof:

If $\lambda>0$ is small enough and Assumption 2 is fulfilled, then formulas
(9) define the resolvent of $\hat{A}$ as an operator that maps the vector
$(\bar{\Delta},\bar{\delta})^{T}\in H$ to ${(\Delta,\delta)^{T}\in H}$.

Let us prove its compactness. For this purpose we consider the Hilbert–Schmidt
norm of the operator ${(a_{ij})=Res(\hat{A})}$,

$\|a_{ij}\|=\sqrt{\sum\limits_{i=1}^{\infty}\sum\limits_{j=1}^{\infty}(a_{ij})^{2}}.$

After performing several estimates based on the Cauchy–Schwarz inequality, we
conclude that $\|Res(\hat{A})\|^{2}$ does not exceed the following value:

$2\sum\limits_{j,i=1}^{\infty}\frac{1}{\omega_{j}^{2}}\left(\frac{1}{\omega_{i}^{2}}\left(\sum\limits_{s,p=1}^{r}\gamma_{s}M^{-1}_{sp}c_{sj}c_{pi}\right)^{2}+2\right).$
(10)

As $|M^{-1}_{sp}|\leq\|M^{-1}\|$, the estimate (10) of $\|Res\hat{A}\|^{2}$ is
finite, while $c_{sj}$ and $\omega_{j}$ satisfy Assumption 2 for every
$s=\overline{1,r}$, $j=1,2,\dots$ . As the Hilbert–Schmidt norm of the
resolvent is finite, $Res(\hat{A}):H\to H$ is a compact operator. Because of
the Dafermos–Slemrod (see [22]), the positive trajectories of system (7) are
precompact in $H$ as the resolvent ${Res(\hat{A})=(\hat{A}-\lambda I)^{-1}}$
is compact for small enough $\lambda>0$. ∎

##  V ASYMPTOTIC STABILITY

###### Assumption 3

The only invariant subspace of ${\rm Ker}\,C$ under the action of the
semigroup $\\{e^{tA}\\}_{t\geq 0}$ is the singleton $\\{0\\}$.

###### Assumption 4

For each $j\in{\mathbb{N}}$, there exists an $s=\overline{1,r}$ such that

$c_{sj}\neq 0.$ (11)

###### Theorem 1

Let Assumptions 1–4 be satisfied, and let the components of the operator $F$
be defined by (6). Then the trivial solution $e=0$ of (7) is asymptotically
stable.

###### Proof:

Consider the following positive definite quadratic form on $H$:

${\cal W}(e)=\sum\limits_{j=1}^{\infty}(\Delta_{j}^{2}+\delta_{j}^{2}),$

and calculate its time derivative along the trajectories of system (7):

$\dot{\cal
W}(e)=-2\sum\limits_{j=1}^{\infty}\sum\limits_{i=1}^{\infty}\sum\limits_{s=1}^{r}c_{si}\,\Delta_{i}\left(f_{js}\Delta_{j}+g_{js}\delta_{j}\right).$

If $f_{js}$ and $g_{js}$ are defined by (6), then

$\dot{\cal
W}(e)=-2\sum\limits_{s=1}^{r}\gamma_{s}\left(\sum\limits_{j=1}^{\infty}c_{sj}\Delta_{j}\right)^{2}\leq
0\;\;\text{for all}\;\left(\begin{array}[]{c}\Delta\\\ \delta\\\
\end{array}\right)\in D(A).$

According to Lyapunov’s theorem, the solution $e(t)\equiv 0$ of (7) is stable.

In order to prove its asymptotic stability, consider the set
$S=\\{e(t):\dot{\cal W}(e)\equiv 0\\}$, where $e(t)$ is the solution of (7).
From the equality

$Ce(t)=0$ (12)

and the structure of $S$, we deduce that $S={\rm Ker}\,C$.

Assumption 3 now appears as the condition of LaSalle’s theorem [23] (see also
[24] for the semigroup formulation), according to which the trivial solution
of the error dynamics (7) is asymptotically stable. ∎

##  VI FLEXIBLE BEAM VIBRATIONS

System (1) represents the operator form of a wide class of mathematical models
of controlled flexible structures. In the sequel, we will consider a
particular case of system (1) resulting from the modal analysis of a vibrating
flexible beam with an attached point mass. A complete description of this
plant can be found in [20]. The equation of motion of this system has been
derived by using Hamilton’s principle, which yields the following relation:

$\displaystyle\int_{0}^{l}\left(\rho\ddot{w}\delta
w+\left(EIw^{\prime\prime}-\sum_{i=1}^{k}u_{i}\psi_{i}\right)\delta
w^{\prime\prime}\right)dx$ (13) $\displaystyle+\left.\left(m\ddot{w}+\varkappa
w-u_{0}\right)\delta w\right|_{x=l_{0}}$ $\displaystyle=0,$

which is assumed to hold for each admissible variation $\delta w(x,t)$ of
class $C^{2}\left([0,l]\times[0,\tau]\right)$, $\tau>0$, satisfying the
boundary conditions

$\delta w|_{t=0}=\delta w|_{t=\tau}=\delta w|_{x=0}=\delta w|_{x=l}=0.$

The function $w(x,t)\in C^{2}\left([0,l]\times[0,\tau]\right)$ represent the
transversal deflection of the beam at a point $x$ and time $t$,
$w(0,t)=w(l,t)=0$; $\rho$, $E$, $I$, $m$, $\varkappa$, and $l$ are positive
mechanical parameters, $l_{0}\in(0,l)$ is the point of the rigid body
attachment, piecewise continuous functions $\psi_{i}(x)$ encode the actuators
placement, ${{\rm supp}\,\psi_{i}\cap\\{0,l_{0},l\\}=\emptyset}$,
$i=\overline{1,k}$.

It has been shown in [21] that the parameters $\omega_{j}$, appearing in the
operator $A$ in (1), are determined in terms of eigenvalues of the spectral
problem

$\begin{gathered}\frac{d^{4}}{dx^{4}}W(x)=\omega^{2}\frac{\rho}{EI}W(x),\quad
x\in(0,l)\setminus\\{l_{0}\\},\\\
W(0)=W(l)=0,\;W^{\prime\prime}(0)=W^{\prime\prime}(l)=0,\;W\in C^{2}[0,l],\\\
W^{\prime\prime\prime}(l_{0}-0)-W^{\prime\prime\prime}(l_{0}+0)=\frac{\varkappa-\omega^{2}m}{EI}W(l_{0}).\end{gathered}$
(14)

Equations (14) are derived by separation of variables in the homogeneous part
of (13).

Precisely, the eigenfrequencies ($\omega_{j}^{2}$ in the current notation) can
be calculated as solutions of the following characteristic equation:

$\Delta(\omega)=0,$ (15)

where

$\displaystyle\Delta(\omega)$
$\displaystyle=\frac{m}{4\mu\rho}\left\\{(\cosh\mu(l-2l_{0})-\cosh\mu
l)\sin\mu l\right.$ $\displaystyle\left.+(\cos\mu(l-2l_{0})-\cos\mu l)\sinh\mu
l\right\\}-\frac{\sin\mu l\sinh\mu l}{\mu^{2}}$
$\displaystyle+\frac{\varkappa}{4EI\mu^{5}}\left\\{(\cosh\mu
l-\cosh\mu(l-2l_{0}))\sin\mu l\right.$
$\displaystyle\left.-(\cos\mu(l-2l_{0})+\cos\mu l)\sinh\mu l\right\\}.$

It has been proved in [19] and [21] that $\omega_{j}$ form an increasing
sequence of distinct positive real numbers, so Assumption 1 is fulfilled.
Moreover (see [20]), the eigenvalues $\omega_{j}^{2}$ of (14) satisfy the
following growth condition: let $Q[a,b)$ denote the number of terms of the
sequence $\\{\omega_{j}\\}_{j\in\mathbb{N}}$ in the interval $[a,b)$, then

$\underset{y\to\infty}{\lim\sup}\;\underset{\tilde{y}\to\infty}{\lim\sup}\;\frac{Q[y,y+\tilde{y})}{\tilde{y}}=0.$

###### Lemma 3

Let $\omega_{j}$ be the solutions of equation (15), ${j=1,2,\dots}$, then the
system of exponents $\\{e^{\pm i\omega_{j}t}\\}_{j=1}^{\infty}$ is linearly
independent on $L^{2}[0,\tau]$ for any $\tau>0$.

The assertion of this lemma follows directly from [25, Theorem 1.2.17].

The eigenfunctions $W_{j}(x)$, $j\in{\mathbb{N}}$ of the corresponding
spectral problem form a linearly independent orthogonal system with respect to
the inner product

$\left<W_{i},W_{j}\right>_{\tilde{H}}=\int_{0}^{l}\rho
W_{i}(x)W_{j}(x)\,dx+mW_{i}(l_{0})W_{j}(l_{0}).$

Let us take the eigenfunctions $W_{1}(x),W_{2}(x),\dots$, corresponding to
eigenvalues $\omega_{1}^{2},\omega_{2}^{2},\dots$, and consider the linear
manifold ${\cal S}={\rm span}\\{W_{1}(x),W_{2}(x),\dots\\}$. Take $w$ and
$\delta w$ from ${\cal S}$, namely,
$w(x,t)=\sum\limits_{i=1}^{\infty}q_{i}W_{i}(x)$ and $\delta w(x,t)=W_{j}(x)$,
and substitute them into (13). Here $q_{i}$ are coefficients of the linear
combination ($q_{i}$ depend on $t$). After performing the integration by parts
and taking into account the orthogonality of the eigenfunctions, we obtain the
following infinite system of ordinary differential equations:

$\ddot{q}_{j}+\omega_{j}^{2}q_{j}=\frac{W_{j}(l_{0})}{\|W_{j}\|^{2}_{\tilde{H}}}u_{0}+\sum_{i=1}^{k}\frac{\int_{0}^{l}\psi_{i}(x)W_{j}^{\prime\prime}(x)\,dx}{\|W_{j}\|^{2}_{\tilde{H}}}u_{i},\;j=1,2,\dots\;.$
(16)

The system (16) in this interpretation represents the orthogonal projection of
the beam-body equations of motion onto the infinite-dimensional linear
manifold ${\cal S}$. In this outline, we treat data $c_{sj}$ as parameters of
the output signals provided by sensors located at points $x=l_{s}$ of the
beam,

$c_{1j}=W_{j}(l_{0}),\quad\text{and}\quad
c_{sj}=W_{j}^{\prime\prime}(l_{s-1}),\quad s=\overline{2,r},$

the control $u_{0}$ is treated as the force applied to the rigid body at
$x=l_{0}$, and $u_{1},\dots,u_{k}$ are treated as actions supplied by $k$
piezoelectric actuators,

$b_{j0}=\frac{W_{j}(l_{0})}{\|W_{j}\|^{2}_{\tilde{H}}},\;\;b_{jp}=\frac{\int_{0}^{l}\psi_{p}(x)W_{j}^{\prime\prime}(x)\,dx}{\|W_{j}\|^{2}_{\tilde{H}}},\;p=1,2,\dots,k.$

Then we denote ${\xi_{j}=\omega_{j}q_{j}}$ and ${\eta_{j}=\dot{q}_{j}}$
respectively, $j\in\mathbb{N}$, and assume that
${\xi=(\xi_{1},\xi_{2},\dots)^{T}\in{\ell^{2}}}$,
${\eta=(\eta_{1},\eta_{2},\dots)^{T}\in{\ell^{2}}}$, which leads to the system
of first-order ordinary differential equations written in the abstract form as
(1)–(2).

Now the assertion of Theorem 1 can be adjusted for the considered flexible
system.

###### Theorem 2

Assume that the components $\omega_{1},\omega_{2},\dots$ of the operator $A$
in (1) are the solutions of equation (15), and the operator $F$ in (7) is
defined by (6) with arbitrary positive gain parameters $\gamma_{s}$,
$s=\overline{1,r}$. If Assumption 4 is satisfied, then the trivial solution
$e=0$ of system (7) is asymptotically stable.

###### Proof:

Instead of step-by-step proof, we will just underline the key difference of
this case from Theorem 1. Note that system (7) degrades into

$\dot{e}(t)=Ae(t)$ (17)

on ${\rm Ker}\,C$.

Suppose initial the conditions

$e(0)=e_{0}\in H.$ (18)

As a consequence of Lemma 1, the Cauchy problem (17), (18) is well-posed on
$H$. The general solution of (17), (18) can be written as

$e(t)=\exp(tA)e_{0}.$ (19)

The substitution of (19) into (12) leads to

$C\exp(tA)e_{0}\equiv 0.$

or,

$\sum\limits_{j=1}^{\infty}c_{sj}\left(\xi_{j}(0)\cos\omega_{j}t+\eta_{j}(0)\sin\omega_{j}t\right)=0,\quad
s=\overline{1,r}.$ (20)

As mentioned in Lemma 3, the system of functions
$\\{\cos\omega_{j}t,\sin\omega_{j}t\\}_{j=1}^{\infty}$ is linearly independent
on $L^{2}[0,\tau]$ for any $\tau>0$, since $\omega_{j}$ satisfy (15). Taking
into account condition (11), the property (20) holds only for
${\xi_{j}(0)=0}$, ${\eta_{j}(0)=0}$, ${\forall\;j=1,2,\dots}$, i.e. for
$e_{0}\equiv 0$. That is, the only solution of system (7) on ${\rm Ker}\,C$ is
the trivial one. So, the only invariant subset of ${\rm Ker}\,C$ under the
action of the semigroup $\\{e^{tA}\\}_{t\geq 0}$ is the point $e=0$.

Thus, we see that Assumption 3 is guaranteed by the asymptotic distribution of
the eigenvalues. So, the conditions of Theorem 1 are satisfied that leads us
to the conclusion about the asymptotic stability of the trivial equilibrium of
(7). ∎

##  VII OBSERVATION ERROR CONVERGENCE

In order to illustrate the observation error dynamics, we perform the
integration of truncated systems (1)–(2), (4), and (7) with the coordinate
indices $j=1,2,...,N$, where $N$ is a given natural number. Let us denote by
$e_{N}(t)$ the solution of the corresponding finite-dimensional system (7)
with $j=1,\dots,N$, and its Euclidean norm by

$\|e_{N}(t)\|=\sqrt{\sum\limits_{j=1}^{N}(\Delta_{j}^{2}+\delta_{j}^{2})}.$

We present the results of numerical simulations for the flexible beam model of
length ${l=1.875\,\text{m}}$ with the rigid body attached at the point
${l_{0}=1.378\,\text{m}}$. This choice of mechanical parameters is described
in [21]. The outputs are assumed to be provided by four piezoelectric sensors
located at ${l_{1}=0.075\,\text{m}}$, ${l_{2}=0.716\,\text{m}}$,
${l_{3}=1.128\,\text{m}}$, and ${l_{4}=1.555\,\text{m}}$. The numerical
integration has been carried out in Maple for the truncated error dynamics
with ${N=6}$, ${N=16}$, and ${N=40}$ modes of vibration. The initial
conditions are taken as
${\Delta_{j}(0)=\delta_{j}(0)=\frac{1}{j\,\omega_{j}}}$.

Figure 1 depicts the first six modal deflections $\Delta_{i}(t)$ with the gain
parameters $\gamma_{i}=6$, $i=1,\dots,6$. Figures 2 and 3 depict,
respectively, $\|e_{16}(t)\|^{2}$ and $\|e_{40}(t)\|^{2}$ with the observer
gain parameters $\gamma_{i}=0.8$, $\gamma_{i}=6$, and $\gamma_{i}=12$,
${i=0,\dots,4}$.

Figure 1: The time plots of $\Delta_{1}(t),\dots,\Delta_{6}(t)$ Figure 2: The
time plot of $\|e_{16}(t)\|^{2}$ Figure 3: The time plot of
$\|e_{40}(t)\|^{2}$

These simulation results illustrate the exponential convergence of the finite-
dimensional error dynamics. We also observe that solutions of the system of
larger dimension ($N=40$) decay slower than the solutions with $N=16$ for
large values of $t$.

##  VIII CONCLUSIONS

We have considered infinite-dimensional Hamiltonian control system (1)–(2) and
derived the Luenberger-type observer in the form (4) that allows to
asymptotically reconstruct the full state of the original system by using a
finite number of outputs. The observer design problem has been solved relying
on Lyapunov’s direct method and the invariance principle, where the energy-
induced Lyapunov functional is taken to measure the observation error. The
proposed observer takes into account the input action and thus admits
arbitrarily large inputs. It also permits unbounded outputs as functions of
time. In the considered flexible beam example, the parameters required for the
observer design are efficiently evaluated in terms of eigenvalues and
eigenfunctions of the corresponding spectral problem.

In this work, we do not raise the question about the decay rate of the
observation error $e(t)$ as $t\to+\infty$. It appears to be exponential for
any finite-dimensional projection of the proposed observation scheme (as
discussed in Section VII), while in the infinite-dimensional case it may not
be exponential (we expect it to be polynomial). The estimation of the
observation error decay rate for the considered class of infinite-dimensional
systems is considered as a topic for future investigation.

## References

  * [1] Luenberger, D. G. (1964). Observing the state of a linear system. IEEE Transactions on Military Electronics, 8(2), 74–80.
  * [2] Zheng, G., Bejarano, F. J., Perruquetti, W., Richard, J. P. (2015). Unknown input observer for linear time-delay systems. Automatica, 61, 35–43.
  * [3] Ahmed-Ali, T., Giri, F., Krstic, M., Kahelras, M. (2018). PDE based observer design for nonlinear systems with large output delay. Systems & Control Letters, 113, 1–8.
  * [4] Yang, K. Y., Wang, J. M. (2019). Pointwise feedback stabilization of an Euler-Bernoulli beam in observations with time delay. ESAIM: Control, Optimisation and Calculus of Variations, 25, 4.
  * [5] Delfour, M. C., Mitter, S. K. (1972). Controllability and observability for infinite-dimensional systems. SIAM Journal on Control, 10(2), 329–333.
  * [6] Kratz, W. (1995). Characterization of strong observability and construction of an observer. Linear Algebra and its Applications, 221, 31–40.
  * [7] Kratz, W., Liebscher, D. (1998). A local characterization of observability. Linear Algebra and its Applications, 269(1-3), 115–137.
  * [8] Demetriou, M. A., Ito, K. (1996). Adaptive observers for a class of infinite dimensional systems. IFAC Proceedings Volumes, 29(1), 5346–5350.
  * [9] Deguenon, A. J., Sallet, G., Xu, C. Z. (2006). A Kalman observer for infinite-dimensional skew-symmetric systems with application to an elastic beam. In Proc. of the Second International Symposium on Communications, Control and Signal Processing.
  * [10] Zuyev, A., Sawodny, O. (2007). Stabilization and observability of a rotating Timoshenko beam model. Mathematical Problems in Engineering, 2007.
  * [11] Liu, Z., Liu, J., Wang, L. (2018). Disturbance observer based attitude control for flexible spacecraft with input magnitude and rate constraints. Aerospace Science and Technology, 72, 486–492.
  * [12] Palomba, I., Richiedei, D., Trevisani, A. (2018). Reduced-order observers for nonlinear state estimation in flexible multibody systems. Shock and Vibration.
  * [13] Mattioni, A., Toledo, J., Le Gorrec, Y. (2020). Observer Based Nonlinear Control of a Rotating Flexible Beam. IFAC-PapersOnLine, 53(2), 7479–7484.
  * [14] Daltin, D. C., Bueno, Á. M., Balthazar, J. M., Serni, P. J. A., Barbosa, R. (2021). State Observer Applied to Position and Vibration Control Using Flexible Link Manipulator. In Vibration Engineering and Technology of Machinery, 281–297. Springer, Cham.
  * [15] Perrusquia, A. (2022). Solution of the linear quadratic regulator problem of black box linear systems using reinforcement learning. Information Sciences, 595, 364–377.
  * [16] Jacob, B., Zwart, H. (2001). Exact observability of diagonal systems with a finite-dimensional output operator. Systems & Control Letters, 43(2), 101–109.
  * [17] Jacob, B., Zwart, H. (2001). Exact observability of diagonal systems with a one-dimensional output operator. International Journal of Applied Mathematics and Computer Science, 11(6), 1277–1283.
  * [18] Russell, D. L., Weiss, G. (1994). A general necessary condition for exact observability. SIAM Journal on Control and Optimization, 32(1), 1–23.
  * [19] A. Zuyev, J. Kalosha (2021). Observer design for a flexible structure with distributed and point sensors. _Proceedings of the Institute of Applied Mathematics and Mechanics of NAS of Ukraine, Vol. 35_(2), 125-136. doi:10.37069/1683-4720-2021-35-9.
  * [20] Kalosha J. I., Zuyev A. L. (2021). Asymptotic stabilization of a flexible beam with an attached mass. _Ukrainian Mathematical Journal, 73_ , 1537–1550.
  * [21] Kalosha J., Zuyev A. and Benner P. (2021). On the Eigenvalue Distribution for a Beam with Attached Masses. _Stabilization of Distributed Parameter Systems: Design Methods and Applications._ _Springer International Publishing,_ 43–56.
  * [22] Dafermos, C. M., Slemrod, M. (1973). Asymptotic behavior of nonlinear contraction semigroups. Journal of Functional Analysis, 13(1), 97–106.
  * [23] LaSalle J. P. (1976). Stability theory and invariance principles. Dynamical systems, 1, 211–222.
  * [24] Zuyev, A. (2003). Partial asymptotic stability and stabilization of nonlinear abstract differential equations. In 42nd IEEE International Conference on Decision and Control (IEEE CDC 2003), 1321–1326.
  * [25] W. Krabs, On Moment Theory and Controllability of One-dimensional Vibrating Systems and Heating Processes, Springer-Verlag, Berlin, Heidelberg, 1992.


# Harmonization Benchmarking Tool for Neuroimaging Datasets

###### Abstract

A major data pre-processing step for large, multi-site studies is to handle
site effects by harmonizing data, generating a dataset that enables more
powerful analyses and more robust algorithms. There is a wide variety of data
harmonization techniques, but there are few tools that streamline the process
of harmonizing data, comparing across techniques, and benchmarking new
techniques. In this paper, we introduce HArmonization BEnchmarking Tool
(HABET), an open source tool for generating harmonized images and evaluating
the performance of different harmonization algorithms. To demonstrate the
capabilities of HABET, we harmonize diffusion MRI images from the Adolescent
Brain and Cognitive Development (ABCD) study using two different approaches,
and we compare their performance.

Index Terms—  Harmonization, Open-source, ComBat, ABCD Study, Diffusion MRI

##  1 Introduction

With the increasing collaboration and pooling of data among neuroimaging
research centers, the ability to reliably conduct mega-analyses is becoming
more important. Traditionally, meta-analyses have often been the only option,
being possible even when only aggregated data is available from each study.
Mega-analysis, in which unaggregated data is pooled across studies, can be
more powerful for analyzing structural neuroimaging data under typical
conditions [1], and it is increasingly an option with multi-site neuroimaging
studies becoming more common [2, 3, 4]. A major challenge when working with
multi-site neuroimaging data is to deal with effects attributable to study
site, such as scanner hardware, reconstruction algorithms, acquisition
protocols, and image quality. These can all affect the final image [5] as well
as downstream results [6, 7]. _Data harmonization_ is the process of removing
these effects. When working to discover biomarkers or construct predictive
models out of multi-site medical image data, harmonization is a critical
preprocessing step to ensure that results can be generalized and applied
across a variety of clinical and research settings.

There is a wide variety of approaches to harmonization. A systematic review of
harmonization literature was conducted in [8], where the authors classified
existing techniques into four broad categories: distribution based techniques,
image processing, synthesis, and invariant feature learning. Distribution
based techniques such as ComBat [9] model the bias between different groups
and subtract that bias out from the original data. Image processing based
techniques such as resampling, image filtering, and normalization, employ
image processing algorithms to harmonize the data. Synthesis and invariant
feature learning methods typically involve deep learning. Synthesis techniques
aim to generate samples from missing groups, similar to the task of style
transfer. For example, this could transfer all of the images from different
sites to match the image data distribution of a particular site. Lastly,
invariant feature learning methods involve encoding image data in some site-
invariant form, and using that encoding instead of the original images in
downstream tasks.

A closely related task to harmonization is _domain adaptation_. This term is
sometimes used to refer to the invariant feature learning [10, 11] or
synthesis [12] categories described above, but it typically refers to the
broader problem of transforming data to match a given data distribution,
beyond the context of multi-center studies.

The Adolescent Brain and Cognitive Development (ABCD) Study [4] is a large,
multi-site, longitudinal study monitoring the brain and cognitive development
of a national cohort of children. Several studies have performed harmonization
on the ABCD dataset. In [13], site effects in the ABCD fMRI images were
examined by training a multinomial logistic classifier to predict scanner from
the imaging data before and after ComBat harmonization [9]. In [14], ABCD 2.0
T1 and T2 weighted images were harmonized with a deep learning approach. In
[15], ABCD DMRI images were harmonized with rotation invariant spherical
harmonics (RISH).

Although many harmonization techniques have been developed, there is a need
for tools that streamline the process of evaluating and comparing the various
techniques on a given dataset. In this work, we introduce the open source
HArmonization BEnchmarking Tool (HABET), and we apply it to investigate site
effects in the ABCD Study dataset. This tool will be available to the research
community as open source software; code can be found at
https://github.com/KitwareMedical/habet. Currently HABET supports ComBat and
global scaling out of the box, and it features a simple framework for adding
new harmonization techniques. In principle, HABET could support techniques
from any of the four aforementioned categories except for invariant feature
learning, as these techniques skip the reconstruction step altogether. More
techniques will be made available in the future. Report generation features
allow researchers to quickly summarize and visualize the performance of
different techniques. Unlike domain adaptation libraries such as [11], our
goal is to provide _evaluation metrics_ for the effectiveness of site effect
removal. A similar work to ours is [12], where the extent of harmonization is
evaluated by comparing features or images from a domain-adapted source dataset
with those from a target dataset. HABET instead focuses its evaluation on
statistical model based detection of site effects with different types of
harmonization, and it accommodates harmonization models that operate on more
than two sites. Our contributions are as follows:

  1. 1.

Introduce the open source framework HABET for harmonization benchmarking.

  2. 2.

Exemplify the capabilities of HABET by analyzing site effects in the ABCD
data, before and after applying two different harmonization techniques.

##  2 Method

HABET provides a simple interface for both harmonizing imaging data and
evaluating the performance of the harmonization technique used. HABET features
two main components, harmonization and report generation; see an overview of
the workflow in Figure 1. These two components are independent, so it is
possible to use HABET only for the purpose of harmonizing images, or for
comparing between different harmonization methods with already harmonized
data.

Both components operate on a given set of images, a mask selecting which part
of the image to process, and a table relating each image to some variables.
This includes the categorical variables whose effect on the image is
undesirable (e.g. site, scanner manufacturer), which we will refer to
collectively as “site.” The table may also include covariates that are of
interest for later analysis and that could be confounded with site. As
previously stated, global scaling and ComBat are supported out of the box,
with more harmonization methods to be added in the future. HABET also features
a simple API for registering custom harmonization techniques, making it
convenient for benchmarking newly developed approaches.

The report generation component aids in quantifying site effects for a given
set of images. This can be executed on a set of images before and after
harmonization, giving the user insight into the performance of a given set of
harmonization methods. At any given voxel location, if there were no site
effects then we would expect the within-site means of the voxel values to be
roughly equal, with some variation due to sampling and other variables. This
hypothesis can be tested via one-way ANOVA: we compute an $F$-statistic and an
associated $p$-value at each voxel, with site as the independent variable and
voxel intensity as the dependent variable. If, for a given voxel, the results
of the ANOVA are found to be significant, then we conclude that there are site
effects at that image location. We may wish to know, for those voxels that
show significant $F$-scores, _which_ sites cause the greatest discrepancy, and
for this we run a $t$-test for each pair of sites. We report, as another
measure of the extent of site effects, the fraction of site pairs that showed
significant effects (similar to the approach in [16]). Bonferroni correction
is applied to the $F$-test $p$-values to account for the number of voxels, and
also to the $t$-test $p$-values to account for the number of voxels and site
pairs. Finally, we want a measure of site effect _magnitude_ when a location
shows significant site effects. For this we report eta-squared values to
indicate the proportion of the variance of each voxel value that can be
attributed to site. Putting it together, each report consists of:

  1. 1.

Site effect significance: a table of voxel-wise $F$-test results from the
ANOVA and an image indicating the subset of the mask where site effects were
detected.

  2. 2.

Site effect size: a table and an image of eta-squared values at each voxel.

  3. 3.

Pairwise site effect significance: a table of voxel-wise $t$-test results and
an image depicting the fraction of significant $t$-tests at each location that
showed significant site effects.

  4. 4.

Pairwise site effect size: a table of Hedges’ $g$ values for voxels that
showed significant site effects.

  5. 5.

A summary table reporting the the total numbers and fractions of $F$-tests and
$t$-tests that were significant.

Fig. 1: The workflow to evaluate a harmonization technique. The input data are
an image dataset (a), an optional mask (b) selecting voxels for inclusion, and
a table (c) associating images to groups (such as study site) over which to
remove effects. Our tool provides two modules: harmonization (d) and report
generation (e). A harmonization technique is configured (d) and an analysis of
site effects is conducted (e), producing reports of the significance and
magnitude of site effects before (f) and after (g) harmonization.

##  3 Results

The ABCD study has 21 participating study sites across the country. A wide
range of data are collected, including structural and functional neuroimaging.
We focus on diffusion tensor imaging (DTI), both as an illustration of our
tool and as a precursor for our future development of tools for white matter
microstructure analysis. Ninety-six diffusion directions were used to acquire
the DMRI images in the ABCD dataset, and the diffusion imaging protocols over
the three different scanner manufacturers were chosen to be as similar as
possible across sites [4]. For details on DMRI acquisition parameters and
preprocessing, see [4]. To examine site effects in the ABCD dataset, we took a
sample of 836 subjects from the baseline year. We ensured that the values of
site, baseline age, sex, scanner hardware, and scanner software version that
occur in the sample were representative of the full ABCD dataset. Deformable
registration to align local structures was achieved with a deep learning
approach based on inverse consistency [17]. All FA images were nonlinearly
registered to the “most typical” subject image, in the sense of [18]: the
subject image selected for use as a template is the one with the minimum mean
displacement when registered to all other images. This approach lets us use a
detailed and population-specific subject template, which improves the quality
of the information available to the registration algorithm and can result in
better alignment. To limit our analysis to white matter, a tissue type
segmentation was computed from the T1-weighted MRI image of the subject image
that was chosen as a common registration target. This was done using FSL for
brain extraction [19] and tissue segmentation [20]. The mask was then eroded
to improve generalizability across the other images.

The images were harmonized using two approaches: ComBat and global scaling.
Reports were generated for each technique and the original images using HABET.

Global scaling fits a simple linear model to the site-averaged images, without
allowing for spatial heterogeneity of site effects. Let $\bar{y{}}{}_{iv}$
denote the value at voxel $v$ averaged over all the images in site $i$, and
let $\bar{y{}}{}_{v}$ denote the value at voxel $v$ averaged over _all_ the
images. The global scaling model, presented in [21], is

$\bar{y{}}{}_{iv}=\theta^{\text{loc}}_{i}+\theta^{\text{scl}}_{i}\,\bar{y{}}_{v}+\epsilon_{iv},$

where $\epsilon_{iv}\sim\mathcal{N}(0,\sigma^{2})$. After fitting the data to
obtain the parameters $\theta^{\text{loc}}_{i}$, $\theta^{\text{scl}}_{i}$,
and $\sigma^{2}$ at each voxel, harmonization proceeds by applying the
transformation
$y\mapsto\frac{y-\hat{\theta}^{\text{loc}}_{i}}{\hat{\theta}^{\text{scl}}_{i}}$
to the values of each image from site $i$. In the terminology of [5], this is
a _global_ as opposed to _voxel-wise_ approach to harmonization. We use global
scaling as a baseline global method to compare with voxel-wise methods like
ComBat.

ComBat, introduced for genomics in [9] and used for imaging data in [21], is
based on the model

$y{}_{ijv}=\alpha_{v}+X_{ij}\beta_{v}+\gamma_{iv}+\delta_{iv}\epsilon_{ijv},$
(1)

where $j$ indexes the subject images. Here the site effects are modeled by the
parameters $\gamma_{iv}$ and $\delta_{iv}$; see [21] for details.

When applying ComBat to the ABCD data, we included age and sex as covariates
in the model (i.e. the matrix $X_{ij}$ in (1)), protecting these variables
from having their effects washed out by the harmonization. We used a
significance level of $\alpha=0.05$ for the ANOVA and $t$-tests.

Table 1 shows summary results for harmonized and unharmonized data: the number
$n_{F}{}$ and fraction $f_{F}{}$ of significant ANOVA tests and the number
$n_{t}{}$ of significant $t$-tests. Figure 2 shows a slice of each output
image from the generated reports. For the unharmonized and global scaling
images, all of the $\eta^{2}$ values from the significant ANOVA tests were
found to be above 0.06. Figure 3 shows the full distribution of $\eta^{2}$
values, regardless of significance, for each harmonization method and for the
unharmonized data.

Harmonization Method | $n_{F}{}$ | $f_{F}{}$ | $n_{t}{}$
---|---|---|---
None (Original Data) | 5123 | 0.43 | 24216
Global Scaling | 3266 | 0.27 | 4887
ComBat | 0 | 0 | N/A
Table 1:  The values of $n_{F}{}$, $f_{F}{}$, and $n_{t}{}$ reported for the
ABCD sample before and after harmonization with global scaling and ComBat.
Note that after harmonization with ComBat, we can no longer detect any site
effects. N/A is reported for the value of $n_{t}{}$ for ComBat because no
$t$-tests were run.

Fig. 2: Axial slices of the spatial site effect report for (a) unharmonized
data and (b) data harmonized wtih global scaling. The values shown in color
are (c) binary ANOVA significance with Bonferroni corrected $\alpha=0.05$, (d)
eta-squared, and (e) the fraction of $t$-tests that were significant with
Bonferroni corrected $\alpha=0.05$. No results are shown here for ComBat
harmonized data, because no voxels were flagged for significant site effects
after ComBat harmonization. The white outline encloses the voxels that were
selected for analysis, and the background image is the T1-weighted MRI
depicting the common space to which all images were registered.

Fig. 3: The estimated distributions of $\eta^{2}$ values for the unharmonized
images and the images harmonized with global scaling and ComBat. Each density
function was estimated from $\eta^{2}$ observations at each of the 11880 voxel
locations under consideration. Note that all voxels are included, rather than
only the ones that had a significant ANOVA.

##  4 Discussion

We have demonstrated how HABET can be used to measure site effects in the ABCD
diffusion MRI data and to compare the performance of different harmonization
techniques.

In the unharmonized images, site effects were detected in nearly half (43%) of
the voxels chosen for analysis. All of the corresponding effect-sizes were
above 0.06, which is typically considered to be a medium effect size. In
short, site effects appear to be common and of considerable magnitude in the
unharmonized images.

Global scaling decreased the percentage of voxels at which site effects were
detected to 27%. This is an improvement from the unharmonized images, but it
is still high. Again, all of the corresponding effect-sizes were above 0.06.
With ComBat, no site effects were detected. In [21], a similar vanishing of
detected site effects was observed with ComBat when analyzing 156 ROI of the
brain.

There are several areas in which our analysis could be improved. First,
Bonferroni correction is a very conservative method of multiple testing
correction, partly because it assumes that the statistical tests at various
voxels are all independent. Image values are typically highly correlated in
space, so the independence assumption is very far from reality. This biases
our analysis heavily towards false negatives (i.e. failing to detect existing
site effects). When we report “0 voxels with significant site effects” after
ComBat harmonization, this could be a statement about the specificity of our
site effect detector as much as it could be a statement about the
effectiveness of ComBat. We can get around this limitation by incorporating
spatially informed techniques for multiple testing correction.

Second, our one-way ANOVA relies on the assumption that the intensity
distribution at a given voxel is normally distributed within each site, with
the same variance across sites. The departure of the data from these
assumptions is currently not reported.

The functionality of HABET may be augmented by introducing better methods for
multiple testing correction, supporting more image formats, and supporting
more harmonization methods out of the box. Finally, we hope to make the tool
more convenient by providing easy access to a default dataset for benchmarking
new techniques.

##  5 Conclusion

In this paper, we presented a new harmonization benchmarking tool, HABET,
which can be used to harmonize a multi-site dataset, compare different
harmonization techniques, and evaluate new harmonization techniques. We
demonstrated the use of HABET by investigating site effects in DMRI images
from the ABCD study dataset and analyzing the effectiveness of ComBat compared
to a simpler harmonization approach.

##  6 Compliance with ethical standards

This research study was conducted retrospectively using data made available by
the NIMH Data Archive. Ethical approval was not required because analyses
using de-identified ABCD Study data are not considered human subjects research
by NIH criteria.

##  7 Acknowledgments

This research was supported in part by the National Institutes of Health under
Award Number R42MH118845. The content is solely the responsibility of the
authors and does not necessarily represent the official views of the National
Institutes of Health.

Data used in the preparation of this article were obtained from the ABCD Study
(https://abcdstudy.org), held in the NIMH Data Archive (NDA). The data used
for the analysis in this work and the raw results obtained can be found at DOI
10.15154/1528197.

## References

  * [1] P. Kochunov et al.,  “Multi-site study of additive genetic effects on fractional anisotropy of cerebral white matter: Comparing meta and megaanalytical approaches for data pooling,”  Neuroimage, vol. 95, pp. 136–150, Jul 2014.
  * [2] F. Alfaro-Almagro et al.,  “Image processing and Quality Control for the first 10,000 brain imaging datasets from UK Biobank,”  Neuroimage, vol. 166, pp. 400–424, 02 2018.
  * [3] B. T. Wyman et al.,  “Standardization of analysis sets for reporting results from ADNI MRI data,”  Alzheimers Dement, vol. 9, no. 3, pp. 332–337, May 2013.
  * [4] Donald J Hagler et al.,  “Image processing and analysis methods for the Adolescent Brain Cognitive Development Study,”  NeuroImage, vol. 202, pp. 116091, Nov. 2019.
  * [5] Maíra Siqueira Pinto et al.,  “Harmonization of Brain Diffusion MRI: Concepts and Methods,”  Frontiers in Neuroscience, vol. 14, pp. 396, May 2020.
  * [6] H. Takao, N. Hayashi, and K. Ohtomo,  “Effect of scanner in longitudinal studies of brain volume changes,”  J Magn Reson Imaging, vol. 34, no. 2, pp. 438–444, Aug 2011.
  * [7] T. Zhu et al.,  “Quantification of accuracy and precision of multi-center DTI measurements: a diffusion phantom and human brain study,”  Neuroimage, vol. 56, no. 3, pp. 1398–1411, Jun 2011.
  * [8] Y. Nan et al.,  “Data harmonisation for information fusion in digital healthcare: A state-of-the-art systematic review, meta-analysis and future research directions,”  Inf Fusion, vol. 82, pp. 99–122, Jun 2022.
  * [9] W. Evan Johnson, Cheng Li, and Ariel Rabinovic,  “Adjusting batch effects in microarray expression data using empirical Bayes methods,”  Biostatistics (Oxford, England), vol. 8, no. 1, pp. 118–127, Jan. 2007.
  * [10] M. Bento, I. Fantini, J. Park, L. Rittner, and R. Frayne,  “Deep Learning in Large and Multi-Site Structural Brain MR Imaging Datasets,”  Front Neuroinform, vol. 15, pp. 805669, 2021.
  * [11] Antoine de Mathelin, François Deheeger, Guillaume Richard, Mathilde Mougeot, and Nicolas Vayatis,  “Adapt : Awesome domain adaptation python toolbox,” 2021.
  * [12] Hao Guan and Mingxia Liu,  “Domainatm: Domain adaptation toolbox for medical data analysis,” 2022\.
  * [13] Dylan M. Nielson et al.,  “Detecting and harmonizing scanner differences in the ABCD study - annual release 1.0,” May 2018,  Pages: 309260 Section: New Results.
  * [14] Siyuan Liu and Pew-Thian Yap,  “Learning Multi-Site Harmonization of Magnetic Resonance Images Without Traveling Human Phantoms,” Sept. 2021,  arXiv:2110.00041 [cs, eess].
  * [15] Suheyla Cetin Karayumak, Fan Zhang, Lauren O’Donnell, and Yogesh Rathi,  “Harmonization of Multi-Site Diffusion Magnetic Resonance Imaging Data From the Adolescent Brain Cognitive Development Study,”  Biological Psychiatry, vol. 91, pp. S84, May 2022.
  * [16] A. Zavaliangos-Petropulu et al.,  “Diffusion MRI Indices and Their Relation to Cognitive Impairment in Brain Aging: The Updated Multi-protocol Approach in ADNI3,”  Front Neuroinform, vol. 13, pp. 2, 2019.
  * [17] H. Greer, R. Kwitt, F. X. Vialard, and M. Niethammer,  “ICON: Learning Regular Maps Through Inverse Consistency,”  Proc IEEE Int Conf Comput Vis, vol. 2021, pp. 3376–3385, Oct 2021\.
  * [18] Stephen M. Smith et al.,  “Tract-based spatial statistics: Voxelwise analysis of multi-subject diffusion data,”  NeuroImage, vol. 31, no. 4, pp. 1487–1505, July 2006.
  * [19] S. M. Smith,  “Fast robust automated brain extraction,”  Hum Brain Mapp, vol. 17, no. 3, pp. 143–155, Nov 2002.
  * [20] Y. Zhang, M. Brady, and S. Smith,  “Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm,”  IEEE Trans Med Imaging, vol. 20, no. 1, pp. 45–57, Jan 2001.
  * [21] Jean-Philippe Fortin et al.,  “Harmonization of multi-site diffusion tensor imaging data,”  NeuroImage, vol. 161, pp. 149–170, Nov. 2017.