2021

[1]Yuya Ohmichi

1]Aviation Technology Directorate, Japan Aerospace Exploration Agency,
Jindaiji-Higashi, Chofu, 182-8522, Tokyo, Japan

2] AdvanceSoft Corporation, Kanda-Surugadai, Chiyoda, 101-0062, Tokyo, Japan

# Time-series image denoising of pressure-sensitive paint data by projected
multivariate singular spectrum analysis

<EMAIL_ADDRESS>Kohmi Takahashi  Kazuyuki Nakakita  [  [

###### Abstract

Time-series data, such as unsteady pressure-sensitive paint (PSP) measurement
data, may contain a significant amount of random noise. Thus, in this study,
we investigated a noise-reduction method that combines multivariate singular
spectrum analysis (MSSA) with low-dimensional data representation. MSSA is a
state-space reconstruction technique that utilizes time-delay embedding, and
the low-dimensional representation is achieved by projecting data onto the
singular value decomposition (SVD) basis. The noise-reduction performance of
the proposed method for unsteady PSP data, i.e., the projected MSSA, is
compared with that of the truncated SVD method, one of the most employed
noise-reduction methods. The result shows that the projected MSSA exhibits
better performance in reducing random noise than the truncated SVD method.
Additionally, in contrast to that of the truncated SVD method, the performance
of the projected MSSA is less sensitive to the truncation rank. The projected
MSSA achieves denoising effectively by extracting smooth trajectories in a
state space from noisy input data. Expectedly, the projected MSSA will be
effective for reducing random noise in not only PSP measurement data, but also
various high-dimensional time-series data.

Graphical Abstract

##  1 Introduction

Understanding the nature of physical phenomena via high-dimensional time-
series images is a fundamental approach employed in scientific research. In
the field of fluid dynamics, the spatiotemporal distribution of pressure and
velocity fields in a flow field can be measured using methods such as
pressure-sensitive paint (PSP) (Liu et al, 2021) and particle-image
velocimetry (Westerweel et al, 2013), for example. However, since experimental
data usually contain unwanted signals, i.e., noise, it is necessary to
properly separate the noise from the signals representing the physical
phenomena.

The PSP measurement, which is the subject of this study, is an optical
measurement technique that has been applied to a variety of flow fields, from
low to hypersonic speeds (Nakakita, 2011; Peng and Liu, 2019; Running et al,
2019), because of its ability to measure pressure with high temporal and
spatial resolution. However, when the pressure fluctuation of the observation
target is small, the PSP image is significantly affected by shot noise (Jiao
et al, 2019; Varigonda and Narayanaswamy, 2021; Liu et al, 2022). Shot noise
is noise that is independent of each pixel, and it generates random
statistical fluctuations that follow the Poisson distribution. Generally,
increasing the light intensity is one method for reducing shot noise. However,
it may be impossible to obtain signals with sufficient intensity from the
unsteady PSP measurement (called fast PSP) because the exposure time is
insufficient. A common technique for reducing shot noise is spatial averaging
(Panda, 2017; Tang et al, 2021). However, this technique does not sufficiently
reduce the noise in some cases, and it reduces the spatial resolution.

Recently, noise-reduction methods based on data-analysis techniques, such as
modal decomposition and sparse sensing, have been explored. The singular value
decomposition (SVD), also called proper orthogonal decomposition (POD)
(Lumley, 1967; Berkooz et al, 1993) or principal component analysis (PCA), are
often used to PSP images to extract coherent patterns or reduce noise
(Pastuhoff et al, 2013; Gordeyev et al, 2014; Peng et al, 2016; Sugioka et al,
2019; Zhou et al, 2022). The SVD is employed with a reduced-rank approximation
(truncated SVD) and noise is removed by reconstructing the PSP images using
only modes that represent signal components. One of the difficulties with the
truncated SVD is extracting only the signal components, which is achieved by
setting the appropriate truncation rank. The choice of the appropriate
truncation rank is nontrivial, and several criteria based on the cumulative
sum of singular values (Pastuhoff et al, 2013), the signal periodicity (Peng
et al, 2016), or reference pressure signals (Sugioka et al, 2019) have been
used. Recently, Gavish and Donoho (2014) proposed the optimal truncation rank
based on asymptotic analysis. However, the appropriate truncation rank is
problem-dependent, and it is difficult to determine the truncation rank
optimally based on a single criterion. In addition, these methods do not
modify the SVD temporal coefficients that represent the temporal behavior of
the modes, even though the coefficients are also affected by noise. Modifying
the temporal coefficients is expected to achieve further noise reduction and
suppress the truncation rank dependence of noise-reduction performance. Wen et
al (2018) proposed a data fusion algorithm that optimizes the temporal
coefficients using a limited number of scattered clean data. Noda et al (2018)
employed the coherent output power between PSP and sound-level meter signals
to extract pressure-distribution patterns relevant to the target periodic
acoustic phenomena. These methods can make effective use of the temporal
behavior information of target phenomena to reduce noise but cannot be applied
to PSP data alone. Recently, Inoue et al (2021) proposed a method that
combines SVD with sparse sensing that can modify the SVD temporal coefficients
using PSP data alone.

Dynamic mode decomposition (DMD) has also been actively studied as a data-
analysis method similar to SVD (Schmid, 2010; Taira et al, 2017). With the DMD
method, one can simultaneously obtain spatial and temporal patterns by
assuming a linear time evolution. Ali et al (2016) applied DMD to PSP images
and extracted pressure-fluctuation patterns as DMD modes. To apply the DMD to
noisy data, several algorithms that reduce the effect of noise have been
proposed (Hemati et al, 2017; Nonomura et al, 2019; Ohmichi et al, 2022).
Ohmichi et al (2022) applied the DMD algorithm for noisy data to fast PSP
images of transonic flow over a swept wing (Ohmichi et al, 2018; Sugioka et
al, 2021), and successfully extracted the pressure distribution related to
transonic buffet phenomenon.

In fluid dynamics research, the time evolution of SVD temporal coefficients is
often analyzed to understand the nature of various phenomena. The phase plot
of the coefficients, which represents the relationship between different SVD
modes, shows a smooth trajectory (Sieber et al, 2016; Edgington-Mitchell et
al, 2018). These researches inspired us to explore a noise-processing method
that combines the low-dimensional representation of data on modal bases with
multivariate singular spectrum analysis (MSSA), which is a multivariate
extension of basic SSA (Broomhead and King, 1986; Hassani and Thomakos, 2010).
Basic SSA is employed for the state-space reconstruction of time-series
signals using time-delay embedding, which is often employed for noise
processing (Golyandina, 2010). By combining MSSA and data representation using
the temporal coefficients and extracting smooth trajectories in the state
space that represents the temporal behavior of physical phenomena, we expect
to achieve effective noise reduction for high-dimensional time-series data.
Notably, a similar approach has been proposed by Chhatkuli et al (2015) in
medical-imaging research. They conducted a time-series forecasting for medical
imaging using the MSSA combined with PCA coefficients, which are identical to
the temporal coefficients of SVD. In this study, we utilized a similar
approach to reduce random noise in high-dimensional time-series data.

The objective of this study is to clarify the effectiveness of a noise-
reduction method that combines MSSA with low-dimensional data representation
via projection onto the SVD basis. For brevity, this method is referred to as
the projected MSSA (PMSSA) in this paper. The projected MSSA is applied to two
test cases. One is the pressure field data of the Kármán vortex street behind
a square cylinder obtained by a numerical fluid simulation, and the other is
the practical unsteady PSP measurement data relevant to a tone trailing-edge
(TE) noise of a two-dimensional airfoil. The noise-reduction performance of
the projected MSSA is compared with that of the truncated SVD (TSVD), which is
one of the most common noise-reduction methods. Additionally, the effects of
the truncation rank used in the projected MSSA and the truncated SVD on the
noise-reduction performance are discussed.

##  2 Approach

###  2.1 Input dataset

The input dataset comprises high-dimensional time-series data, such as time-
series spatial-distribution data. We assume that each instantaneous data,
denoted as a snapshot in this paper, has $d$ elements ($d$ is the number of
pixels of the images) and $j$-th snapshot $\bm{x}(j)$ can be represented as
$\bm{x}(j)=[x_{1,j},\cdots,x_{d,j}]^{\mathsf{T}}$, where $j=1,\cdots,m$ and
$m$ is the number of snapshots. Thus, the input dataset can be represented as
the matrix $X$ whose columns are $\bm{x}(j)$,

${X}=\begin{bmatrix}\bm{x}(1)&\cdots&\bm{x}(m)\end{bmatrix}=\begin{bmatrix}x_{1,1}&\cdots&x_{1,m}\\\
\vdots&\ddots&\vdots\\\ x_{d,1}&\cdots&x_{d,m}\end{bmatrix}.$ (1)

That is, $X$ is a $d\times m$ matrix, where the row and column indices
correspond to the spatial and temporal distributions, respectively.

###  2.2 Noise reduction by truncated SVD

The input data matrix $X$ can be factorized by SVD to

$X=USV^{\mathsf{T}}=\sum_{i=1}^{r_{X}}\bm{u}_{i}\sigma_{i}\bm{v}_{i}^{\mathsf{T}},$
(2)

where $U$ and $V$ are the orthogonal matrices, and $S$ is the diagonal matrix
containing the singular values. The left and right singular vectors,
$\bm{u}_{i}$ and $\bm{v}_{i}$, are the $i$-th columns of $U$ and $V$,
respectively; $\sigma_{i}$ is the $i$-th singular value, and $r_{X}$ is the
rank of $X$. The singular values are arranged in descending order. The
$\bm{u}_{i}$ and $\bm{v}_{i}$ represent a spatial mode and its temporal
behavior, respectively. If the input dataset $X$ is extremely large, hindering
the application of the batch SVD algorithm, online algorithms (Ross et al,
2008; Halko et al, 2011; Ohmichi, 2017) can be applied.

In the noise-reduction method by the truncated SVD, we assume that the noise
can be separated as the component corresponding to small singular values
($\sigma_{i}<\sigma_{r}$).

$X=\sum_{i=1}^{r}\bm{u}_{i}\sigma_{i}\bm{v}_{i}^{\mathsf{T}}+\sum_{i=r+1}^{r_{X}}\bm{u}_{i}\sigma_{i}\bm{v}_{i}^{\mathsf{T}}.$
(3)

The first and second terms in the right hand side of Eq. 3 represent the
signal and noise components, respectively. That is, the noise-reduced data
$X_{svd}$ are obtained as

$X_{svd}=U_{r}S_{r}V_{r}^{\mathsf{T}}=\sum_{i=1}^{r}\bm{u}_{i}\sigma_{i}\bm{v}_{i}^{\mathsf{T}}.$
(4)

The truncation rank $r$ is often selected based on the cumulative sum of
singular values or prior information, such as the signal periodicity.

Note that the SVD is identical to the POD method; however, in the POD method,
mean subtraction from the input dataset is essential.

###  2.3 Noise reduction by the projected MSSA

In this study, we examine the noise-reduction performance of the projected
MSSA, which is composed of the MSSA and the low-dimensional representation
using the SVD basis.

#### Step 1: Projecting data onto the subspace

Using SVD, the input data $X$ can be projected onto the subspace spanned by
the first $r$ singular vectors $U_{r}$:

$\tilde{X}=U_{r}^{\mathsf{T}}X=S_{r}V_{r}^{\mathsf{T}}.$ (5)

The $\tilde{X}$ represents the low-dimensional representation of the input
data and is a $r\times m$ matrix. We denote the $i$-th row and element of
$\tilde{X}$ as $\bm{\tilde{x}}_{i}$ and ${\tilde{x}}_{i,j}$, respectively.

$\tilde{X}=\begin{bmatrix}\bm{\tilde{x}}_{1}\\\ \vdots\\\
\bm{\tilde{x}}_{r}\end{bmatrix}=\begin{bmatrix}\tilde{x}_{1,1}&\cdots&\tilde{x}_{1,m}\\\
\vdots&\ddots&\vdots\\\ \tilde{x}_{r,1}&\cdots&\tilde{x}_{r,m}\end{bmatrix}.$
(6)

As shown in Eqs. 4 and 5, $\bm{\tilde{x}}_{i}$ (denoted as the temporal
coefficient hereinafter) represents the temporal change in the expansion
coefficients for $\bm{u}_{i}$. For example, if $r=r_{X}$, the $j$-th snapshot
can be expressed using $\tilde{x}_{i,j}$ as

${\bm{x}}(j)=\tilde{x}_{1,j}\bm{u_{1}}+\tilde{x}_{2,j}\bm{u}_{2}+\cdots+\tilde{x}_{r,j}\bm{u}_{r}.$
(7)

Thereafter, the noise reduction is implemented by applying MSSA to
$\tilde{X}$. In other words, this algorithm applies MSSA to the SVD-based
denoised data via projected low-dimensional representation. The MSSA consists
of the following steps.

#### Step 2: Embedding

The SSA reduces noise by extracting smooth trajectories in a state space. The
state-space reconstruction is performed by time-delay embedding as follows. We
define the trajectory matrix $\tilde{T}_{i}$ of $\bm{\tilde{x}}_{i}$ as

${\tilde{T_{i}}}=\begin{bmatrix}\tilde{x}_{i,1}&\tilde{x}_{i,2}&\cdots&\tilde{x}_{i,K}\\\
\tilde{x}_{i,2}&\tilde{x}_{i,3}&\cdots&\tilde{x}_{i,K+1}\\\
\vdots&\vdots&\ddots&\vdots\\\
\tilde{x}_{i,L}&\tilde{x}_{i,L+1}&\cdots&\tilde{x}_{i,m}\end{bmatrix}.$ (8)

Here, $L$ and $K$ represent the window length and the window width,
respectively, and $K=m-L+1$. Each column of $\tilde{T}_{i}$ is a partial time-
series with the length of $L$ extracted from the time series
$\bm{\tilde{x}}_{i}$ and represents a state-space vector converted from
$\bm{\tilde{x}}_{i}$ by time-delay embedding. $\tilde{T}_{i}$ has the same
elements along the antidiagonals. The multivariate trajectory matrix
${\tilde{T}}$ is constructed by stacking each trajectory matrix
$\tilde{T}_{i}$:

${\tilde{T}}=\begin{bmatrix}\tilde{T}_{1}\\\ \vdots\\\
\tilde{T}_{r}\end{bmatrix}.$ (9)

$\tilde{T}$ is a $rL\times K$ matrix.

#### Step 3: SVD and grouping

In general SSA procedures, the trajectory matrix undergoes SVD, and the
decomposed components are grouped by some criteria. Here, the grouping is
performed following the same procedure employed for noise reduction using the
truncated SVD method (Section 2.2). That is, we decompose the trajectory
matrix $\tilde{T}$ by SVD as ${\tilde{T}}={P}{D}{Q}^{\mathsf{T}}$, and assume
that the noise components can be separated as the components corresponding to
small singular values. Then, the noise-reduced trajectory can be extracted as

$\hat{T}=\begin{bmatrix}\hat{T}_{1}\\\ \vdots\\\
\hat{T}_{r}\end{bmatrix}={P}_{r_{mssa}}{D}_{r_{mssa}}{Q}_{r_{mssa}}^{\mathsf{T}}$
(10)

Here, we set the truncation rank $r_{mssa}$ to $r_{mssa}=r$.

#### Step 4: Diagonal averaging

Since the trajectory matrix has elements corresponding to temporal
coefficients of the same snapshot along the antidiagonals (Eq. 8), the noise-
reduced trajectory matrix $\hat{T}_{i}$ can be transformed into the temporal
coefficient $\bm{\hat{x}}_{i}$ by obtaining the average along the
antidiagonals:

${\hat{x}_{i,j}}=\sum_{(a,b)\in A_{j}}\frac{(\hat{T}_{i})_{a,b}}{n(A_{j})},$
(11)

and

$\hat{X}=\begin{bmatrix}\bm{\hat{x}}_{1}\\\ \vdots\\\
\bm{\hat{x}}_{r}\end{bmatrix}=\begin{bmatrix}\hat{x}_{1,1}&\cdots&\hat{x}_{1,m}\\\
\vdots&\ddots&\vdots\\\ \hat{x}_{r,1}&\cdots&\hat{x}_{r,m}\end{bmatrix}.$ (12)

Here, $(\hat{T}_{i})_{a,b}$ is the element in the $a$-th row and $b$-th column
of $\hat{T}_{i}$. $A_{j}$ represents a set of indices $(a,b)$, comprising the
antidiagonals of $\hat{T}_{i}$, i.e., $A_{j}=\\{(a,b):a+b-1=j,1\leq a\leq
L,1\leq b\leq K\\}$. $n(A_{j})$ is the number of elements in $A_{j}$.

#### Step 5: Reconstructing the noise-reduced data

Finally, the noise-reduced time-series data $X_{mssa}$ can be obtained in a
manner similar to Eq. 7 as

$X_{mssa}=U_{r}\hat{X}.$ (13)

##  3 Numerical experiments using fluid-simulation data

###  3.1 Input data

First, the noise-reduction performances of TSVD and PMSSA were evaluated for
the pressure field of a Kármán vortex street behind a square cylinder obtained
by spatial two-dimensional fluid simulation. The Mach number of the flow field
was 0.2, and the Reynolds number was 120, defined using the length $l$ of the
side of the square cylinder and the velocity $U_{\infty}$ of the uniform flow
as reference values. For details on this fluid simulation, please refer to the
previous study (Ohmichi et al, 2022). The pressure-coefficient distribution of
the wake region, i.e., $x/l=[0,~{}10]$ and $y/l=[-5,~{}5]$, sampled into
equally spaced $101\times 101$ grid points was used as the input dataset. The
center of the cylinder was located at the origin $(x/l,y/l)=(0,0)$. The number
of snapshots was $m=1000$, and the time interval between the snapshots was
$\Delta t=0.125l/U_{\infty}$. The shot noise could be modeled as random normal
noise. To investigate the noise-reduction performance, random normal noise
with a standard deviation $\sigma$ was added to each snapshot as the
observation noise.

###  3.2 Results and discussion

Figure 1 shows the noise-reduced images by TSVD and PMSSA. The top row shows
the input-pressure distribution with noise values of $\sigma=0,~{}0.08$, and
0.64. The middle and bottom rows show the results of the noise reduction for
the input data with $\sigma=0.08$ and $0.64$, respectively. The window length
in PMSSA was set to $L=90$, and the truncation ranks were set to $r=5$ and 11.
Figure 1 indicates that when the noise was small $\sigma=0.08$, the
differences between TSVD and PMSSA and the different $r$ results were visually
small. The result of PMSSA had a slightly smoother distribution at high noise
level $\sigma=0.64$ than that of TSVD.

Figure 1: Input-pressure distribution images with different noise levels on
the top row, and the noise-reduced images with TSVD and PMSSA on the middle
($\sigma=0.08$) and bottom ($\sigma=0.64$) rows.

To quantitatively evaluate the noise-reduction performance, we define the
relative error as follows:

$Relative~{}error=\frac{{\left\lVert
X_{rec}-X_{0}\right\rVert}_{F}}{{\left\lVert X_{0}\right\rVert}_{F}}$ (14)

$X_{rec}$ represents the data denoised by TSVD and PMSSA (i.e., $X_{svd}$ and
$X_{mssa}$), and $X_{0}$ is the original data without noise ($\sigma=0$).

The effect of the window length $L$ on the relative error in PMSSA is shown in
Fig. 2. The appropriate value of $L$ may depend on the nature of the signal
and noise; however, Golyandina (2010) reported that $L=m/2$ is a good choice
in many cases for the basic SSA. Figure 2 shows that $L=m/2$ exhibited a small
error here as well. However, in the interval of
$30~{}(\approx\sqrt{m})<L<500~{}(=m/2)$, the difference in error was
relatively small. The size of the trajectory matrix is $rL\times(m-L+1)$;
thus, the smaller the $L$, the lower the computational cost. An empirical
guideline for achieving the appropriate $L$ involves choosing a value in the
order of $\sqrt{m}$. In this problem, we set $L=90~{}(\approx 3\sqrt{m})$.

Figure 2: Effect of the window length $L$ on the relative error.
$\sigma=0.64.$

Figure 3 shows the change in error for different noise levels at each
truncation rank. For TSVD, the difference in the truncation rank had a
significant impact on the error. The error was the smallest at $r=5$ at high
noise levels. At low noise levels, the errors were smaller at $r=11$ and
$r=21$ than that of $r=5$ because even the signal component was truncated.
PMSSA is less sensitive to the truncation rank, and a larger truncation rank
than that for TSVD can be applied even at high noise levels. When $r=11$ and
21, small errors were achieved at all noise levels. At sufficiently low noise
levels, the errors of TSVD and PMSSA were comparable.

Figure 3: Effect of the truncation rank $r$ on the relative error in TSVD and
PMSSA.

Figure 4 shows the spatial modes $\bm{u}_{i}$ extracted by SVD at each noise
level. The first mode $\bm{u}_{1}$ represents the mean field, while the second
and third modes represent the pressure patterns indicating Kármán vortex
shedding. The fourth/fifth and sixth/seventh mode pairs represent the second
and third harmonics of the Kármán vortex mode, respectively, while the sixth
and seventh modes at $\sigma=0.64$ did not exhibit coherent patterns because
the noise level was extremely high.

Figure 4: Spatial modes $\bm{u}_{i}$ extracted by SVD at different noise
levels.

Next, we examined the reproducibility of the time-series change. Figure 5
shows the frequency distribution for the pressure at the point, ($x/l,~{}y/l$)
= (3.0, 0.5), behind the square cylinder. For each truncation rank, the
frequency distributions for different noise levels are compared. Figure 6
shows the comparison in the time domain. From these figures, we can see that
PMSSA was superior to TSVD in removing the noise components under all
conditions. Figure 5 shows that TSVD could extract periodic signals, although
the noise remained at frequencies different from those of the periodic signal
components. Conversely, PMSSA reduced noise over a wide frequency range while
retaining the periodic signal representing a Kármán vortex shedding
phenomenon. This noise-reduction performance of PMSSA was not significantly
affected by the truncation rank whereas the remaining noise increased with the
truncation rank for TSVD. Figures 5a and 5b show that at low noise levels
($\sigma=0.08$) and with an overly small truncation rank ($r=5$), even the
signal components, which appeared as harmonic peaks at $r=11$, were truncated.
At high noise levels, this harmonic peak was not captured even at $r=11$
because the amplitude of the harmonic signal was smaller than the noise floor.
This is a limitation of SVD-based methods: since both PMSSA and TSVD
reconstruct the signal using SVD modes, if the SVD does not capture a target
phenomenon, these methods cannot reproduce the phenomenon.

Figure 5: Comparison of the pressure signals at the point, ($x/l,~{}y/l$) =
(3.0, 0.5), behind the square cylinder in the frequency domain.  Figure 6:
Comparison of the pressure signals at the point, ($x/l,~{}y/l$) = (3.0, 0.5),
behind the square cylinder in the time domain.

Finally, we compared the temporal coefficients $\bm{\tilde{x}}_{i}$ and
$\bm{\hat{x}}_{i}$ (Eqs. 6 and 12). As explained in Section 2.3, these are the
results of the noise reduction by TSVD and PMSSA, expressed as the expansion
coefficient for $\bm{u}_{i}$. Figure 7a shows the time history of the
coefficient for $\bm{u}_{1}$, which is the mean mode. Evidently, the
coefficient for the mean mode fluctuated significantly in the TSVD, while the
fluctuation was suppressed in the PMSSA. Figures 7b and 7c show the phase
plots for the coefficients of $\bm{u}_{2},~{}\bm{u}_{3}$ and
$\bm{u}_{4},~{}\bm{u}_{5}$, respectively. These pairs represent the Kármán
vortex and its second-harmonic oscillation, respectively. The figures clearly
show that while PMSSA was able to extract smooth trajectories representing the
relationship between the modes, the TSVD result retained a significant amount
of noise in the trajectories. These results indicate that PMSSA was superior
in extracting the smooth trajectories of physical signals from noise-
contaminated signals.

Figure 7: Comparison of the temporal coefficients $\bm{\tilde{x}}_{i}$ and
$\bm{\hat{x}}_{i}$ of TSVD and PMSSA noise-reduced data for fluid-simulation
data analyses. $\sigma=0.64$.

##  4 Demonstration using PSP data

###  4.1 Input data

As a demonstration using practical experimental data, we applied the PMSSA and
TSVD to PSP data containing large noise. The analyzed data were the time-
series measurement data of the wing surface pressure distribution for the tone
TE noise phenomenon on the NACA0012 airfoil (Nakakita, 2011, 2013). The
pressure distribution was measured using anodized aluminum PSP, which is one
of the fast PSPs. The low-turbulence calibration wind tunnel at JAXA Chofu
Aerospace Center was employed for this experiment. Figure 8 shows the
experimental model and the measured pressure-distribution image. The main flow
velocity was 28.0 m/s, the angle of attack was $-1.5^{\circ}$, and the
pressure measurement area was a circular region on the pressure side of the
wing (indicated by an arrow in Fig. 8a). The sampling rate of the pressure-
distribution image was 10 kHz, and the signals below 300 Hz were removed by a
frequency filter to eliminate the dominant low-frequency noise, which was not
random noise. A $2\times 2$ spatial averaging was performed, and the number of
pixels after the spatial averaging was $d=23917$. The number of snapshots was
$m=10000$. Figure 8b shows the pressure-distribution image after applying the
frequency filter and spatial averaging. Since the shot noise was large
compared to the physical signal, coherent patterns could not be observed
visually. However, Nakakita (2011, 2013) succeeded in extracting a coherent
pressure pattern at approximately 918 Hz relevant to the tone TE noise by
cross-spectrum analysis from this data (shown in Fig. 8c). In the present
analysis, we apply PMSSA to this data to demonstrate its effectiveness in
reducing shot noise. For details on this experiment, please refer to Nakakita
(2011, 2013).

Figure 8: (a) NACA0012 airfoil model and measurement area for the PSP
measurement. (b) Instantaneous pressure-distribution image contaminated by
shot noise. (c) Coherence map of the 918 Hz signal taken by Nakakita (2013).

###  4.2 Results and discussion

First, we consider the results of the mode decomposition by SVD. Figure 9
shows the spatial mode $\bm{u}_{i}$ obtained by SVD on the upper row and the
frequency distribution of the temporal coefficient $\bm{\tilde{x}}_{i}$
corresponding to each mode on the lower row. The most dominant mode,
$\bm{u}_{1}$, had a noisy spatial distribution with a slightly coherent
structure. The frequency with the largest amplitude was 916 Hz, which
corresponds to the fundamental frequency of the TE noise; however, the
amplitude was similar in all frequency ranges, indicating that this mode
behaved almost randomly in time. Conversely, the second and third modes
($\bm{u}_{2}$ and $\bm{u}_{3}$) exhibited a clear wavy pattern, and the
frequency distribution indicated that they displayed periodic oscillations at
916 Hz. The fourth mode ($\bm{u}_{4}$) exhibited a similar wavy pattern with a
peak frequency of 916 Hz, although different spatial patterns from the second
and third modes were present. In the frequency distribution of this mode, the
amplitude of $O(100)$ Hz was slightly higher, and the different patterns shown
in the spatial distribution could be considered to have $O(100)$ Hz
fluctuations. Although not shown in the figure, the fifth and subsequent modes
did not have a clear peak frequency.

Figure 9: Spatial modes $\bm{u}_{i}$ extracted by SVD (upper row), and the
frequency distributions of their temporal coefficients $\bm{\tilde{x}}_{i}$
(lower row).

Based on the above discussion of the SVD mode, the truncation ranks for TSVD
and PMSSA were set to $r=4$ to capture the 916 Hz fluctuation pattern. In
PMSSA, the results at $r=11$ were also compared to examine the effect of the
truncation rank. The window length was set to $L=3\sqrt{m}=300$.

Figure 10 shows the pressure-distribution image after the noise removal by the
TSVD and PMSSA methods. In the image before denoising (Fig. 8b), shot noise
appeared in high levels, and coherent patterns were not observed. Contrarily,
after denoising by TSVD and PMSSA, in the image shown in Fig. 10, wavy
patterns corresponding to the periodic fluctuation of 916 Hz appeared.
However, in the TSVD results, the pressure values of the entire region
fluctuated significantly, which is attributed to the random-noise variation of
$\bm{u}_{1}$ shown in Fig. 9. In PMSSA, only the temporal change in the wavy
pattern was clearly extracted. It was also observed that PMSSA was less
sensitive to the truncation rank $r$, as it successfully removed the noise
even when $r=11$, which included an excessive number of modes to extract the
916 Hz oscillation.

Figure 10: Reconstructed (i.e., denoised) pressure-distribution images
obtained by TSVD ($r=4$), PMSSA ($r=4$), and PMSSA ($r=11$).

In Fig. 11, the comparison of the pressure fluctuations at the position
indicated by the arrow in Fig. 8b is shown. The result of the pressure
measurement by the Kulite pressure transducer is also shown for comparison.
From the frequency distribution shown in Fig. 11a, it is evident that the peak
amplitudes at 916 Hz for both methods of TSVD and PMSSA were comparable to
that by the pressure transducer. However, for the TSVD, the noise levels in
other frequency ranges were high. PMSSA could significantly reduce the noise
in both $r=4$ and 11. From the comparison of the pressure fluctuations in the
time domain in Fig. 11b, it was confirmed that the denoised signal with PMSSA
exhibited less noise than that with TSVD, and a smooth temporal signal was
extracted. As shown in Fig. 11a, the second-harmonic peak frequency of
approximately 1840 Hz was captured by the pressure transducer, but not by TSVD
and PMSSA. This is because the amplitude of the pressure fluctuation at 1840
Hz was smaller than that of the shot noise generated in the PSP data, and the
mode decomposition by SVD could not extract this phenomenon.

Figure 11: Comparisons of the denoised pressure signals and pressure
transducer signal in the (a) frequency and (b) time domain.

Figure 12 shows a comparison of the temporal coefficients
($\bm{\tilde{x}}_{i}$ and $\bm{\hat{x}}_{i}$) obtained by the TSVD and PMSSA
methods. The time variation of the temporal coefficient for the first mode
($\bm{u}_{1}$) in Fig. 12a clearly shows that the PMSSA reduced the large
random fluctuations remaining in the TSVD result. In addition, in Fig. 12b,
the phase plot represented by the coefficients for the second and third modes
($\bm{u}_{2}$ and $\bm{u}_{3}$) clearly shows that the PMSSA extracted a
smoother trajectory than TSVD. In other words, PMSSA can identify the modes
behaving randomly and modes representing physical phenomena via modifying
their trajectories. The effective denoising is achieved by reconstructing PSP
data based on the identified physically meaningful SVD modes and their
modified trajectories.

Figure 12: Comparison of the temporal coefficients $\bm{\tilde{x}}_{i}$ and
$\bm{\hat{x}}_{i}$ of the TSVD and PMSSA noise-reduced data for the PSP data
analysis.

##  5 Conclusion

A numerical method for reducing random noise from high-dimensional time-series
data was investigated for denoising unsteady PSP data. Unsteady PSP data are
significantly affected by random noise generated by shot noise. Here, we
demonstrated the effectiveness of a noise-reduction method, i.e., the
projected MSSA. The projected MSSA comprises a combination of MSSA and low-
dimensional data representation. For the low-dimensional representation, we
projected the data onto the subspace spanned by the SVD basis. The pressure
distributions for Kármán vortex shedding obtained by numerical fluid
simulations and PSP data on the NACA0012 airfoil relevant to the tone TE noise
were analyzed. The results clearly showed several advantages of the projected
MSSA. The projected MSSA exhibited better noise-reduction performance than the
conventional truncated SVD method in both analyses. In addition, the projected
MSSA was less sensitive to the truncation rank than the truncated SVD and
exhibited better performance in a wide range of truncation rank values, since
the projected MSSA can remove random variations from the SVD temporal
coefficients. Furthermore, the projected MSSA does not require prior
information such as noise amplitude or pressure tap data as reference signals
to identify noise components. These advantages are achieved by extracting
smooth trajectories of temporal coefficients in a state space by using time-
delay embedding. It should be noted, however, that since the projected MSSA
uses the spatial modes from SVD for low-dimensional projection, signals not
captured by SVD modes cannot be recovered. The projected MSSA can be combined
with mode decomposition methods other than SVD, and the development of mode
decomposition methods with better detection limits is desired.

## Acknowledgments

The authors thank Hiroki Iwamoto for his technical support regarding the
processing of the PSP data.

## Declarations

### Ethical Approval

Not applicable.

### Competing interests

The authors have no competing interests to disclose.

### Authors’ contributions

YO contributed to conceptualization, data curation, funding acquisition,
investigation, methodology, software, validation, visualization, writing –
original draft, and writing – review and editing; KT contributed to
investigation, software, validation, writing – original draft, and writing –
review and editing; KN contributed to investigation, resources, and writing –
review and editing.

### Funding

This work was supported in part by the Japan Society for the Promotion of
Science (JSPS) KAKENHI (grant no. 20K14958).

### Availability of data and materials

The data that support the findings of this study are not publicly available.
Data are however available from the authors upon reasonable request and with
permission of Japan Aerospace Exploration Agency.

## References

  *   * Ali et al (2016) Ali MY, Pandey A, Gregory JW (2016) Dynamic Mode Decomposition of Fast Pressure Sensitive Paint Data. Sensors (Basel, Switzerland) 16(6):862. 10.3390/s16060862
  * Berkooz et al (1993) Berkooz G, Holmes P, Lumley JL (1993) The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows. Annual Review of Fluid Mechanics 25(1):539–575. 10.1146/annurev.fl.25.010193.002543
  * Broomhead and King (1986) Broomhead DS, King GP (1986) Extracting qualitative dynamics from experimental data. Physica D: Nonlinear Phenomena 20(2):217–236. 10.1016/0167-2789(86)90031-X
  * Chhatkuli et al (2015) Chhatkuli RB, Demachi K, Miyamoto N, et al (2015) Dynamic Image Prediction Using Principal Component and Multi-Channel Singular Spectral Analysis: A Feasibility Study. Open Journal of Medical Imaging 05(03):133–142. 10.4236/ojmi.2015.53017
  * Edgington-Mitchell et al (2018) Edgington-Mitchell DM, Weightman JL, Honnery DR, et al (2018) Sound production by shock leakage in supersonic jet screech. In: 2018 AIAA/CEAS Aeroacoustics Conference. AIAA 2018-3147, 10.2514/6.2018-3147
  * Gavish and Donoho (2014) Gavish M, Donoho DL (2014) The optimal hard threshold for singular values is $4/\sqrt{3}$. IEEE Transactions on Information Theory 60(8):5040–5053. 10.1109/TIT.2014.2323359
  * Golyandina (2010) Golyandina N (2010) On the choice of parameters in singular spectrum analysis and related subspace-based methods. Statistics and Its Interface 3(3):259–279. 10.4310/SII.2010.v3.n3.a2
  * Gordeyev et al (2014) Gordeyev S, De Lucca N, Jumper EJ, et al (2014) Comparison of unsteady pressure fields on turrets with different surface features using pressure-sensitive paint. Experiments in Fluids 55(1):1661. 10.1007/s00348-013-1661-9
  * Halko et al (2011) Halko N, Martinsson PG, Tropp JA (2011) Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions. SIAM Reviews 53:217–288. 10.1137/090771806
  * Hassani and Thomakos (2010) Hassani H, Thomakos D (2010) A review on singular spectrum analysis for economic and financial time series. Statistics and Its Interface 3:377–397. 10.4310/SII.2010.v3.n3.a11
  * Hemati et al (2017) Hemati MS, Rowley CW, Deem EA, et al (2017) De-biasing the Dynamic Mode Decomposition for Applied Koopman Spectral Analysis of Noisy Datasets. Theoretical and Computational Fluid Dynamics 31(4):349–368. 10.1007/s00162-017-0432-2
  * Inoue et al (2021) Inoue T, Matsuda Y, Ikami T, et al (2021) Data-driven approach for noise reduction in pressure-sensitive paint data based on modal expansion and time-series data at optimally placed points. Physics of Fluids 33(7):077,105. 10.1063/5.0049071
  * Jiao et al (2019) Jiao L, Chen Y, Wen X, et al (2019) Resolving vortex-induced pressure fluctuations on a cylinder in rotor wake using fast-responding pressure-sensitive paint. Physics of Fluids 31(5):055,106. 10.1063/1.5092944
  * Liu et al (2021) Liu T, Sullivan J, Asai K, et al (2021) Pressure and Temperature Sensitive Paints, 2nd edn. Springer International Publishing, Cham, Switzerland, 10.1007/978-3-030-68056-5
  * Liu et al (2022) Liu X, Qin C, Tang Y, et al (2022) Resolving dynamic features of kilohertz pressure fluctuations using fast-responding pressure-sensitive paint: Measurement of inclined jet impingement. Experiments in Fluids 63(4):72. 10.1007/s00348-022-03419-4
  * Lumley (1967) Lumley JL (1967) The Structure of Inhomogeneous Turbulent Flows. In: Yaglom AM, Tatarski VI (eds) Atmospheric turbulence and radio propagation. Nauka, Moscow, p 166–178
  * Nakakita (2011) Nakakita K (2011) Unsteady Pressure Measurement on NACA0012 Model Using Global Low-Speed Unsteady PSP Technique. In: 41st AIAA Fluid Dynamics Conference and Exhibit. AIAA 2011-3901, 10.2514/6.2011-3901
  * Nakakita (2013) Nakakita K (2013) Phase Delay and Correlation Area Detection of Unsteady Pressure Field Using Unsteady PSP Measurement. In: AIAA Ground Testing Conference. AIAA 2013-3124, 10.2514/6.2013-3124
  * Noda et al (2018) Noda T, Nakakita K, Wakahara M, et al (2018) Detection of small-amplitude periodic surface pressure fluctuation by pressure-sensitive paint measurements using frequency-domain methods. Experiments in Fluids 59(6):94. 10.1007/s00348-018-2550-z
  * Nonomura et al (2019) Nonomura T, Shibata H, Takaki R (2019) Extended-Kalman-filter-based dynamic mode decomposition for simultaneous system identification and denoising. PLOS ONE 14(2):e0209,836. 10.1371/journal.pone.0209836
  * Ohmichi (2017) Ohmichi Y (2017) Preconditioned dynamic mode decomposition and mode selection algorithms for large datasets using incremental proper orthogonal decomposition. AIP Advances 7(7):075,318. 10.1063/1.4996024, https://arxiv.org/abs/1704.03181
  * Ohmichi et al (2018) Ohmichi Y, Ishida T, Hashimoto A (2018) Modal Decomposition Analysis of Three-Dimensional Transonic Buffet Phenomenon on a Swept Wing. AIAA Journal 56(10):3938–3950. 10.2514/1.J056855
  * Ohmichi et al (2022) Ohmichi Y, Sugioka Y, Nakakita K (2022) Stable dynamic mode decomposition algorithm for noisy pressure-sensitive-paint measurement data. AIAA Journal 60(3):1–6. 10.2514/1.J061086
  * Panda (2017) Panda J (2017) Experimental Verification of Buffet Calculation Procedure Using Unsteady Pressure-Sensitive Paint. Journal of Aircraft 54(5):1791–1801. 10.2514/1.C033917
  * Pastuhoff et al (2013) Pastuhoff M, Yorita D, Asai K, et al (2013) Enhancing the signal-to-noise ratio of pressure sensitive paint data by singular value decomposition. Measurement of Science and Technology 24(7):075,301. 10.1088/0957-0233/24/7/075301
  * Peng and Liu (2019) Peng D, Liu Y (2019) Fast pressure-sensitive paint for understanding complex flows: From regular to harsh environments. Experiments in Fluids 61(1):8. 10.1007/s00348-019-2839-6
  * Peng et al (2016) Peng D, Wang S, Liu Y (2016) Fast PSP measurements of wall-pressure fluctuation in low-speed flows: Improvements using proper orthogonal decomposition. Experiments in Fluids 57(4):45. 10.1007/s00348-016-2130-z
  * Ross et al (2008) Ross DA, Lim J, Lin RS, et al (2008) Incremental Learning for Robust Visual Tracking. International Journal of Computer Vision 77(1-3):125–141. 10.1007/s11263-007-0075-7
  * Running et al (2019) Running CL, Sakaue H, Juliano TJ (2019) Hypersonic boundary-layer separation detection with pressure-sensitive paint for a cone at high angle of attack. Experiments in Fluids 60(1):23. 10.1007/s00348-018-2665-2
  * Schmid (2010) Schmid PJ (2010) Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanicsurnal of Fluid Mechanics 656:5–28. 10.1017/S0022112010001217
  * Sieber et al (2016) Sieber M, Paschereit CO, Oberleithner K (2016) Spectral proper orthogonal decomposition. Journal of Fluid Mechanics 792:798–828. 10.1017/jfm.2016.103
  * Sugioka et al (2019) Sugioka Y, Hiura K, Chen L, et al (2019) Unsteady pressure-sensitive-paint (PSP) measurement in low-speed flow: Characteristic mode decomposition and noise floor analysis. Experiments in Fluids 60(7):108. 10.1007/s00348-019-2755-9
  * Sugioka et al (2021) Sugioka Y, Nakakita K, Koike S, et al (2021) Characteristic unsteady pressure field on a civil aircraft wing related to the onset of transonic buffet. Experiments in Fluids 62(1):20. 10.1007/s00348-020-03118-y
  * Taira et al (2017) Taira K, Brunton SL, Dawson STM, et al (2017) Modal Analysis of Fluid Flows: An Overview. AIAA Journal 55(12):4013–4041. 10.2514/1.J056060
  * Tang et al (2021) Tang L, Hand L, Murakami D, et al (2021) Unsteady Pressure-Sensitive-Paint Shot Noise Reduction. In: AIAA AVIATION 2021 FORUM. AIAA 2021-2579, 10.2514/6.2021-2579
  * Varigonda and Narayanaswamy (2021) Varigonda SV, Narayanaswamy V (2021) Methodology to image the panel surface pressure power spectra in weakly coupled fluid/structure interactions. Experiments in Fluids 62(11):237. 10.1007/s00348-021-03314-4
  * Wen et al (2018) Wen X, Liu Y, Li Z, et al (2018) Data mining of a clean signal from highly noisy data based on compressed data fusion: A fast-responding pressure-sensitive paint application. Physics of Fluids 30(9):097,103. 10.1063/1.5046681
  * Westerweel et al (2013) Westerweel J, Elsinga GE, Adrian RJ (2013) Particle image velocimetry for complex and turbulent flows. Annual Review of Fluid Mechanics 45(1):409–436. 10.1146/annurev-fluid-120710-101204
  * Zhou et al (2022) Zhou L, Wang S, Song J, et al (2022) Study of internal time-resolved flow dynamics of a subsonic fluidic oscillator using fast pressure sensitive paint. Experiments in Fluids 63(1):17. 10.1007/s00348-021-03370-w


# Data Diversity based Metamorphic Relation Prioritization

Madhusudan Srinivasan  Montana State University
Bozeman, United States
<EMAIL_ADDRESS>Upulee Kanewala  University of North Texas
Florida, United States
<EMAIL_ADDRESS>
###### Abstract

An oracle is a mechanism to decide whether the outputs of the program for the
executed test cases are correct. For machine learning programs, such oracle is
not available or too difficult to apply. Metamorphic testing is a testing
approach that uses metamorphic relations, which are necessary properties of
the software under test to help verify the correctness of a program.
Prioritization of metamorphic relations helps to reduce the cost of
metamorphic testing [1]. However, prioritizing metamorphic relations based on
code coverage is often not effective for prioritizing MRs for machine learning
programs, since the decision logic of a machine learning model is learned from
training data, and 100% code coverage can be easily achieved with a single
test input. To this end, in this work, we propose a cost-effective approach
based on diversity in the source and follow-up data set to prioritize
metamorphic relations for machine learning programs. We show that the proposed
data diversity-based prioritization approach increase the fault detection
effectiveness by up to 40% when compared to the code coverage-based approach
and reduce the time taken to detect a fault by 29% when compared to random
execution of MRs. Overall, our approach leads to saving time and cost during
testing.

###### Index Terms:

Metamorphic Testing, Metamorphic Relations, Machine Learning

##  I Introduction and Motivation

Machine learning (ML) is increasingly deployed in large-scale software and
safety-critical systems due to recent advancements in deep learning and
reinforcement learning. The software applications powered by ML are applied
and utilized in our daily lives; from finance, and energy, to health and
transportation. AI-related accidents are already making headlines, from
inaccurate facial recognition systems causing false arrests to unexpected
racial and gender discrimination by machine learning software [2] [3] [4].
Also, a recent incident with an Uber autonomous car resulted in the death of a
pedestrian [5]. Thus, ensuring the reliability and conducting rigorous testing
of machine learning systems is very vital to prevent any future disasters.

One of the critical components of software testing is the test oracle. The
mechanism to determine the correctness of test case output is known as test
oracle [6]. A program can be considered a non-testable program if there does
not exist a test oracle or too difficult to determine whether the test output
is correct [7]. Hence, those programs suffer from a problem known as the test
oracle problem. ML-based systems exhibit non-deterministic behavior and reason
in probabilistic manner. Also, the output is learned and predicted by a
machine learning model as result do not have an expected output to compare
against the actual output. As a result, correctness of the output produced by
machine learning-based applications cannot be easily determined and suffer
from test oracle problem [8].

Metamorphic testing (MT) is a property-based software testing technique that
is used to alleviate the test oracle problem. A central component of MT is a
set of Metamorphic Relations (MRs), which are necessary properties of the
target function or algorithm in relation to multiple inputs and their expected
outputs [9] [10]. In recent years especially, MT has attracted a lot of
attention and detected a large number of real-life faults in various domains.
For example, Murphy et al [11] identified six metamorphic relations and real
bugs were identified on three machine learning tools. In a related project,
Xie et al [12], [13] used metamorphic testing for the detection of faults in
K-nearest neighbours and Naive Bayes classifiers. The results revealed that
both algorithms violated some of the metamorphic relations. Also, some real
faults were detected in the open–source machine learning tool Weka [14]. Zhang
et al [15] applied MT to test DNN based autonomous driving systems. MT was
able to find inconsistent behaviors under different road scenes for real-world
autonomous driving models.

Previous work [16] has shown that several MRs detect the same type of fault
and differ only in fault detection effectiveness. Also, each MR in a machine
learning application can have multiple source and follow-up test cases, as a
result, the execution time of the MRs increases drastically. The source and
follow-up test cases in a MR represented as training data set could take
several days to train and build a machine learning model for a large training
data set. Hence, prioritization of MRs is important for conducting effective
metamorphic testing on these programs.

Prior work [1] introduced approaches to prioritize MRs based on code coverage
and fault-based approach. The fault-based prioritization approach involves
using mutant killed information of each MR to prioritize MR. The code coverage
approach uses statement or branch coverage information for prioritizing MRs.
Further, the code coverage-based approach does not work effectively for
machine learning-based applications since the data used for training
determines the program logic. Machine learning program is usually just a
sequence of library function invocations and as a result 100% code coverage
can be easily achieved when at least one test input is processed. To this end,
in this work, we propose methods for _MR prioritization_ for machine learning
programs. The MR prioritization method proposed in this work qualitatively
identifies the diversity between the source and follow-up training data set in
a MR. Moreover, our results indicate that MR prioritization reduces the number
of MRs needed to execute on SUT and reduces the time taken to detect a fault
up to 29%.

We make the following contributions in this work:

  1. 1.

Propose a novel data diversity based automated method where we proposed four
metrics to prioritize _metamorphic relations_ based on diversity between the
source and follow-up data set in a MR.

  2. 2.

Evaluate the effectiveness of our proposed MR prioritization methods on four
machine learning algorithm implementations and perform a comparison against
two baseline approaches (1) Random baseline: represents the current practice
of executing the MRs in a random order, (2) Code coverage-based ordering:
represents the MR ordering based on the statement or branch coverage
information collected when testing the previous version of the SUT as show in
the prior work [1].

##  II Proposed Approach

Figure 1: Steps for MR prioritization of ML Programs

In this section, we discuss our proposed method for prioritization of MRs.
Figure 1 shows the steps for prioritizing MRs. In this approach, we prioritize
the MRs as follows:

  1. 1.

Let the set of source test cases in a MR used for testing the SUT be the
_prioritizing source test cases_ $(T_{sp})$.

  2. 2.

Let the set of follow-up test cases in a MR used for testing the SUT be the
_prioritizing follow-up test cases_ $(T_{fp})$.

  3. 3.

For $(T_{sp})$ and $(T_{fp})$, apply the proposed metric: rule based, anomaly
detection, clustering based and data distribution. The metrics are discussed
in detail in section II-A.

  4. 4.

For $(T_{sp})$ and $(T_{fp})$, apply the proposed metric: rule based, anomaly
detection, clustering based and data distribution and get the MR diversity
value of each MR based on the formulas provided in section II-A.

  5. 5.

MR diversity value calculated using the metrics for each MR is normalized
using the following steps:

    1. (a)

Let the set of MR diversity value of MRs be _$DMR_{v}$_.

    2. (b)

Identify minimum and maximum value in $DMR_{v}$. Let minimum value be
$Min_{v}$ and maximum value be $Max_{v}$.

    3. (c)

For each MR diversity value $v\in DMR_{v}$, generate the normalized MR
diversity value using the formula below:

$N_{MR}=\frac{v-Min_{v}}{Max_{v}-Min_{v}}$ (1)

Where $N_{MR}$ represent the normalized MR diversity value.

  6. 6.

Prioritize the normalized MR diversity values and select the top _n_ MRs from
the prioritized ordering to execute based on the resources available for
testing.

Our proposed metrics are explained in detail below.

###  II-A Dataset Diversity Approach

In this work, we propose four metrics to prioritize metamorphic relations
based on diversity between source and follow-up data set in a MR. Our
intuition behind the metrics is that greater the diversity between the source
and follow-up data set, then greater the quality of MR and performance of the
ML model.

The proposed metrics are explained in detail below.

####  II-A1 Rule Based Classifier

Rule-based classifiers are just another type of classifier which makes the
class decision by using various if-then rules. We apply a rule based
classifier (CN2) to determine diversity between the source and follow-up data
set in a metamorphic relation. CN2 is an instance of sequential covering
algorithm family based on Rule-Based classification used for learning set of
rules.The basic idea here is to learn one rule, remove the data that it
covers, then repeat the same process. We hypothesize that if greater the
diversity in the classification rules between the source and followup data
set, then greater the performance of the model and MR. We apply the following
steps.

  1. 1.

Let the set of source test cases used for testing SUT be the _prioritizing
source test cases_ $(T_{sp})$.

  2. 2.

Let the set of follow-up test cases of the SUT be the _prioritizing follow-up
test cases_ $(T_{fp})$.

  3. 3.

Apply Sequential Covering Algorithm (CN2) on $T_{sp}$ and $T_{fp}$ and obtain
the rules.The CN2 induction algorithm is a learning algorithm for induction of
simple and comprehensible rules.

  4. 4.

Remove same set of rules covered in both $T_{sp}$ and $T_{fp}$.

  5. 5.

Calculate the total rules identified in $T_{sp}$ and $T_{fp}$. Let
$Total_{sr}$ and $Total_{fr}$ represent the total rules in $T_{sp}$ and
$T_{fp}$ respectively.

  6. 6.

We calculate the diversity of MR using the formula below:

$\begin{split}\text{$R_{MR}$}=\text{$\mid(Total_{sr})-(Total_{fr})\mid$}\end{split}$
(2)

Where $R_{MR}$ represent the diversity of the MR.

####  II-A2 Anomaly Detection

Anomaly detection (or outlier detection) is the identification of rare items,
events or observations which raise suspicions by differing significantly from
the majority of the data. We apply anomaly detection to identify outliers in
the source and followup data set in a metamorphic relation. We hypothesize
that greater the diversity in outliers detected between the source and
followup, greater the possibility of detecting fault in the machine learning
model. We apply the following steps.

  1. 1.

Let the set of source test cases used for testing SUT be the _prioritizing
source test cases_ $(T_{sp})$.

  2. 2.

Let the set of follow-up test cases of the SUT be the _prioritizing follow-up
test cases_ $(T_{fp})$.

  3. 3.

Apply K Nearest Neighbour (KNN) algorithm to $T_{sp}$ and $T_{fp}$ and
identify outlier instances.

  4. 4.

Remove identical outliers detected between $T_{sp}$ and $T_{fp}$.

  5. 5.

Calculate the total outliers in $T_{sp}$ and $T_{fp}$ respectively. Let
$Total_{so}$ and $Total_{fo}$ be the sum of outlier detected in $T_{sp}$ and
$T_{fp}$ respectively.

  6. 6.

We calculate the diversity of MR using the formula below:

$\begin{split}\text{$O_{MR}$}=\text{$\mid(Total_{so})-(Total_{fo})\mid$}\end{split}$
(3)

Where $O_{MR}$ represent the diversity of the MR.

####  II-A3 Clustering-based

In this metric, we partition the input space ie. all possible inputs into
different region and determine how the data is distributed in the space. To
achieve this, we use clustering methods which simply try to group similar
patterns into clusters whose members are more similar to each other based on
some distance measure than to members of other clusters. We apply the
following steps.

  1. 1.

Let the set of source test cases used for testing SUT be the _prioritizing
source test cases_ $(T_{sp})$.

  2. 2.

Let the set of follow-up test cases of the SUT be the _prioritizing follow-up
test cases_ $(T_{fp})$.

  3. 3.

Apply K mean algorithm to $T_{sp}$ and $T_{fp}$ and identify clusters. Let
clusters of $T_{sp}$ and $T_{fp}$ be $S_{cl}$ and $F_{cl}$.

  4. 4.

Find the distance between the clusters by calculating the euclidean distances
between the clusters in $S_{cl}$. Let the total distance between the clusters
be $Ts_{cds}$.

  5. 5.

Determine the size of each cluster in $S_{cl}$. Let the total size of the
clusters in $S_{cl}$ be $Ts_{sc}$.

  6. 6.

Find the distance within the cluster by calculating the average of the
distances from observations to the centroid of each cluster in $S_{cl}$. Let
average distance within the cluster in $S_{cl}$ be $As_{cd}$.

  7. 7.

Apply the steps 4 to 6 to $F_{cl}$ respectively. Let $Tf_{cds}$ be the total
distance between the clusters, $Tf_{sc}$ be the total size of the clusters and
$Af_{cd}$ be the average distance within the cluster.

$\begin{split}\text{$R_{MR}$}=\text{$\mid(Total_{sr})-(Total_{fr})\mid$}\end{split}$
(4)

  8. 8.

We calculate the diversity of MR using the formula below:

$\begin{split}\begin{split}\text{$C_{MR}$}=\text{$\mid(Ts_{cds}+Ts_{sc}+As_{cd})-(Tf_{cds}+Tf_{sc}+Af_{cd})\mid$}\end{split}\end{split}$
(5)

Where $C_{MR}$ is the diversity of MR.

####  II-A4 Data Distribution

Data used in training machine learning model often form very similar patterns.
For example, a lot of data were grouped around the middle values, with fewer
observations at the outside edges of the distribution. These patterns are
known as distributions, because they describe how the data are distributed
across the range of possible values. In this metric, we used shape and spread
of data distribution of dataset to find diversity of metamorphic relations.
The spread of distribution involves calculating range, variance and standard
deviation of the source and follow-up data set. Similarly, shape of
distribution involves calculating skewness and kurtosis of the source and
follow-up data set of the MR. We hypothesize that larger the diversity in the
data distribution between the source and follow-up data set of a MR, then
greater the quality of MR. We apply the following steps.

  1. 1.

Let the set of source test cases used for testing SUT be the _prioritizing
source test cases_ $(T_{sp})$.

  2. 2.

Let the set of follow-up test cases of the SUT be the _prioritizing follow-up
test cases_ $(T_{fp})$.

  3. 3.

Calculate the skewness and kurtosis for $T_{sp}$. Let $ST_{sk}$ be the sum of
skewness and kurtosis for $T_{sp}$.

  4. 4.

Calculate range, variance and standard deviation for $T_{sp}$. Let $ST_{rvs}$
be the sum of variance, range and standard deviation for $T_{sp}$.

  5. 5.

Apply step 3 and 4 to $T_{fp}$. Let $FT_{sk}$ and $FT_{rvs}$ be the shape and
spread of distribution for $T_{fp}$.

  6. 6.

We calculate the diversity of MR using the formula below:

$\begin{split}\text{$DD_{MR}$}\begin{split}=\text{$\mid(ST_{sk}+ST_{rvs})-(FT_{sk}+FT_{rvs})\mid$}\end{split}\end{split}$
(6)

Where $DD_{MR}$ is the diversity of MR.

##  III Evaluation

In this work, we plan to evaluate the utility of the developed prioritization
approaches on the following aspects: (1) Fault detection effectiveness of MR
sets, (2) Effective number of MRs required for testing, and (3) Time taken to
detect a fault. We evaluate the effectiveness of the proposed approaches with
(1) a _random baseline_ : this represents the current practice of executing
source and follow-up data set of the MRs in random order. (2) _code-coverage
based ordering_ : this represents the MR order based on the statement and
branch-coverage information of the MRs. The approach selects MR which has
covered highest number of statements or branches and place it in the
prioritized MR ordering. The process continues until all the possible
statement or branches are covered [1].

We conducted experiments to find answers for the following research questions:

  1. 1.

Research Question 1 (RQ1): _Are the proposed MR prioritization approach based
on data diversity more effective than the random baseline?_

  2. 2.

Research Question 1 (RQ2): _Are the proposed MR prioritization approach more
effective than the code coverage-based prioritization?_

###  III-A Evaluation Procedure

In order to answer the above research questions, we carry out the following
validation procedure similar to previous work [1]:

  1. 1.

We generate a set of mutants for the SUT. These mutants will be referred to as
the _validation set of faults_ $(F_{v})$. Then from the generated mutants, we
remove mutants that give exceptions or mutants that do not finish execution
within a reasonable amount of time from further consideration.

  2. 2.

We use the method described in Section II to obtain data diversity-based MR
ordering. Then we applied the obtained MR ordering to $F_{v}$ and log the
mutant killing information.

  3. 3.

Creating the random baseline: We generate 100 random MR orderings and apply
each of those orderings to $F_{v}$. The mutant killing information for each of
those random orderings is logged. The average mutant killing rates of these
100 random orderings were computed to get the fault detection effectiveness of
the random baseline.

  4. 4.

Creating the Code Coverage-based ordering: Apply greedy approach on the
statements or branches covered to prioritize the MRs as follows [1]:

    1. (a)

Select the MR that covers the highest number of statement or branches in the
program and place it in the prioritized MR ordering.

    2. (b)

Remove the statements and branches covered by the MR.

    3. (c)

Repeat the previous two steps until all the statements and branches are
covered.

###  III-B Evaluation Measures

We used the following measures to evaluate the effectiveness of the MR
orderings generated by our proposed metrics:

  1. 1.

To measure the _fault detection effectiveness_ of a given set of MRs, we use
the percentage of mutants killed by those MRs.

  2. 2.

To calculate the _effective MR set size_ we used the following approach:
typically, the fault detection effectiveness of MT increases as the number of
MRs used for testing increases. However, after a certain number of MRs, the
rate of increase in fault detection slows due to factors such as redundancy of
MRs. Therefore when there is no significant increase in the fault detection
between two MR sets of consecutive sizes of size $m$ and $m+1$, where the MR
set of size $m+1$ is created by adding one MR to the MR set set of size $m$,
the effective MR set size can be determined. That is, if the difference in
fault detection effectiveness of MR set of size $m$ and MR set of size $m+1$
is less than some threshold value, $m$ would be the effective MR set size that
should be used for testing. Determining this threshold value should be done
considering the critical nature of the SUT. In this work, we used two
threshold values of 5% and 2.5% as used in previous related work for
determining the oracle data set size [17].

  3. 3.

We used the following approach to find the _average time taken to detect a
fault_ : for each killable mutant $m$ in $F_{v}$, we computed the time taken
to kill the mutant ($t_{m}$) by computing the time taken to execute the source
and follow-up test cases of the MRs in a given prioritized order until $m$ is
killed (here it is assumed that the source and follow-up test cases for each
MR are executed sequentially). Then the average time taken to detect a fault
is computed using the following formula:

$\frac{\sum t_{m}}{\\#killable\,mutants\,in\,F_{v}}$

  4. 4.

We calculate _Average percentage of fault detected_ (APFD) developed by Elbaum
et al. [18] [19] [20] that measures the average rate of fault detection per
percentage of test suite execution. The APFD is calculated by taking the
weighted average of the number of faults detected during the execution of MRs.
APFD can be calculated using the formula below:

$\begin{split}\text{APFD}=\frac{1-MRF_{1}+MRF_{2}+...+MRF_{m}}{nm}+\frac{1}{2n}\end{split}$
(7)

Where $MRF_{i}$ represents the fault $F_{i}$ detected by the metamorphic
relations under evaluation, m represents the number of faults present in the
SUT and n represents the total number of MRs used.

###  III-C Subject Programs and MRs

We applied the above-mentioned validation procedure on the following four
machine learning programs for evaluating our proposed MR prioritization
methods.

  * •

IBk111http://weka.sourceforge.net/doc.dev/weka/classifiers/lazy/IBk.html:
K-nearest neighbours (KNN) classifier implementation in the Weka machine
learning library [21]. The input to IBk is a training data set, and a test
data set to represent in .arff format. The output is the classification
predictions made on the instances in the test data set.

  * •

Linear
Regression222https://weka.sourceforge.io/doc.dev/weka/classifiers/functions/LinearRegression:
linear regression is a linear approach that models the relationship between
the input variables and the single output variable. The subject program is the
linear regression implementation in weka library. Weka latest stable version
is 3.8.5, which consists of the multiple linear regression, named
LinearRegression() class. The input to linear regression program is a training
data set, and a test data set to represent in .arff format. The output is the
prediction made on the instances in the test data set.

  * •

Naive
Bayes333https://weka.sourceforge.io/doc.dev/weka/classifiers/bayes/NaiveBayes.html:
A Naive Bayes classifier is a probabilistic machine learning model that’s used
for classification task. The subject program is the Naive Bayes algorithm
implementation, named NaiveBayes() class in Weka library. The input to Naive
Bayes program is a training data set, and a test data set to represent in
.arff format. The output is a classification made on the instances in the test
data set.

  * •

Convolutional Neural
Network444https://javadoc.io/doc/org.deeplearning4j/deeplearning4j-nn/latest/index.html:
Convolutional Neural Network (CNN) is a class of Artificial Neural Network and
most commonly applied to analyze visual imagery. The subject program is the
CNN algorithm implementation in Deeplearning4j library. The input to the CNN
is MNIST handwritten digit data set and the output is to classify the image of
a handwritten digit into one of 10 classes representing integer values from 0
to 9.

###  III-D Metamorphic Relations

For conducting MT on IBk and Naive Bayes, we used 11 MRs developed by Xie et
al. [12]. These MRs are developed based on the user expectations of supervised
classifiers. These MRs modify the training and testing data so that the
predictions do not change between the source and follow-up test cases.
Similarly, for testing Linear Regression system, we used 11 MRs developed by
Luu et al. [22]. The properties of the linear regression is related to
addition of data points, the rescaling of inputs, the shifting of variables,
the reordering of data, and the rotation of independent variables. These MRs
were grouped into 6 categories and derived from the properties of the targeted
algorithm. For testing convolutional neural network, we used 11 MRs developed
by Ding et al. [23]. The MRs were developed on three levels: system level,
data set level and data item level. The MRs in the data set level are created
on the relation of the classification accuracy of reorganized training data
sets. The MRs in the data item level are created on the relation of the
classification accuracy of reproduced individual images.

###  III-E Source and Follow-up Test Cases

As we described before, MT involves the process of generating source and
follow-up test cases based on the MRs used for testing.

IBk uses a training data set to train the k-nearest neighbor classifier, and a
test data set is used to evaluate the performance of the trained classifier.
After executing source test case, the output will be the class labels
predicted for the test data set, which is a value from the set {0,1,2,3,4}.
The attribute values are randomly selected within the range [0, 100]. The
values of the class label are randomly selected from the set {0,1,2,3,4,5}.
The size of the training testing data sets ranges within [0, 500]. To generate
the follow-up test cases using these source test cases, we applied the input
transformations described in the MRs.

Linear Regression uses a training data set to train the model, and a test data
set is used to evaluate the performance of the trained model. We obtained the
training data from machine learning repository
555https://archive.ics.uci.edu/ml/datasets/YearPredictionMSD as source data
set. The training data contains 515345 instances and 90 attributes. To
generate the follow-up data set using these source data set, we applied the
input transformations described in the MRs. Similarly for Naive Bayes, we
obtained the training data from machine learning repository
666https://archive.ics.uci.edu/ml/datasets/Adult as source data set. The
training data contains 48842 instances and 14 attributes. To generate the
follow-up data set using these source data set, we applied the input
transformations described in the MRs. For testing convolutional neural
network, we obtained the training data 777http://yann.lecun.com/exdb/mnist/ as
source data set. The training dataset contains 60000 handwritten digits and
10000 images for testing. To generate the follow-up test cases using these
source test cases, we applied the input transformations described in the MRs.

###  III-F Mutant Generation

For each subject program, we aimed at developing generated mutant set to be
used as $F_{v}$. For this, we used three automated mutation tools:
$\mu$Java888https://github.com/jeffoffutt/muJava, PIT999http://pitest.org/ and
Major101010http://mutation-testing.org/. For generating mutants using
$\mu$Java, we used all the method level mutation operators [24]. With PIT
mutation tool [25] and the Major mutation tool [26] we used all the available
mutation operators provided by the tools. All the mutants the were generated
using these tools were _first order mutants_ , where a mutation operator is
used to make a single modification to the source code to generate the mutant.
Table I shows the number of mutants used for evaluation.

.

TABLE I: Mutants generated for version $v_{k}$ of the SUTs Version | Subject | #Major | #PIT | # $\mathbf{\mu}$Java
---|---|---|---|---
3.8.6 | IBk | 1021 | N/A | 621
3.8.6 | Linear Regression | 377 | 73 | N/A
3.8.6 | Naive Bayes | 452 | 452 | N/A
0.9.1 | CNN | 92 | N/A | N/A
TABLE II: Validation Setup Subject | $\mathbf{{T_{sp}}}$ | $\mathbf{{F_{v}}}$
---|---|---
IBk | Test1 | Major
Linear Regression | YearPrediction | Major
Naive Bayes | Adult | PIT
CNN | MNIST | PIT

##  IV Result and Analysis

In this section, we discuss our experimental results and provide answers to
the two research questions that we listed in Section III. For each subject
program, we carried out the validation procedure described in Section III-A
using the setup described in Table II. In this setup, we used the generated
mutant set for evaluating the prioritized MR ordering as $F_{v},$ and the
source test cases as $T_{sp}$. For example, in Table II for IBk, Test1 refers
to the randomly generated data set that we described in Section III-E.
Similarly, for CNN, MNIST refers to the 60000 images used for training and
testing. Adult and YearPrediction refers to the training data set used for
Naive Bayes and Linear Regression program respectively. Figure 2 shows the
average fault detection effectiveness for evaluation runs described above vs.
the MR set size used for testing, for each subject program. We also plot the
percentage of faults detected by the random baseline, statement and branch
coverage based MR ordering for comparison. Results for each research question
is provided below.

(a) IBk

(b) CNN

(c) Naive Bayes

(d) Linear Regression

Figure 2: Fault detection effectiveness for IBk, Linear Regression, Naive
Bayes and Convolutional Neural Network

###  IV-A RQ1: Comparison of the MR Prioritization approaches against the
Random Baseline

####  IV-A1 Fault detection effectiveness

We formulated the following statistical hypothesis to answer RQ1 in the
context of fault detection effectiveness:

  * $H_{1}$:

For a given MR set of size m, the fault detection effectiveness of the MR set
produced by data diversity-based approach is higher than that of the random
baseline.

The null hypothesis $H_{0x}$ for each of the above defined hypothesis $H_{X}$
is that the data-diversity based approaches perform equal to or worse than the
random baseline.

In Table III to VI, we list the relative improvement in average fault
detection effectiveness for the data diversity-based method over the currently
used random approach of prioritizing MRs for Linear Regression, Naive Bayes,
IBk and CNN respectively. To evaluate the above hypotheses, we used the paired
permutation test since it does not make any assumptions about the underlying
distribution of the data [27]. We apply the paired permutation test to each MR
set size for each of the subject programs with $\alpha$= 0.05, and the
relative improvements that are significant are marked with a *.

As shown in Table III, for IBk, rule based metric shows improvements in
average fault detection effectiveness over the random approach for all the MR
set sizes except for the last MR set size. Particularly, among the different
MR set sizes, the increase in fault detection percentage varies from 1.23% to
29.5%. Therefore, we reject the null hypotheses $H_{01}$ for IBk for rule
based, Anomaly, clustering based and data distribution metric. Similarly for
Linear Regression, Table IV indicates that the anomaly based, clustering based
and data distribution based metrics show consistent significant improvements
in fault detection over random. Therefore we can reject $H_{01}$ in general
for Regression. The fault detection varies from 2.5% to 30% for the three
metrics. Table V shows the relative improvement in the average fault detection
percentage for Naive Bayes. Our proposed metrics shows improvements for MR set
sizes $m=1$ to $m=5$. Also overall, the relative improvement across MR
prioritization methods varies between 0% - 40%. Therefore we reject $H_{01}$
in general for Naive Bayes. Table VI shows the relative improvement in the
average fault detection percentage for CNN. Clustering-based metric shows
improvements for all MR set sizes except MR set size $m=5$. Also overall, the
relative improvement across MR prioritization methods varies between 0% -66%.
Therefore we can reject $H_{01}$ in general for CNN.

Therefore our data diversity based approach outperforms the random approach
for all our subject programs.

TABLE III: Relative improvement in fault detection effectiveness of the data diversity-based approach compared to random baseline for IBk MR Set-Size |  Rule-based |  Anomaly-based |  Clustering-based |  Data Distribution
---|---|---|---|---
1 |  29.54%* |  29.54%* |  25%* |  20.45%*
2 |  17.54%* |  8.77%* |  5.26%* |  19.29%*
3 |  11.11%* |  3.17%* |  9.52%* |  9.52%*
4 |  5.88%* |  1.47%* |  4.41%* |  2.94%*
5 |  5.63%* |  4.22%* |  7.04%* |  2.81%*
6 |  5.40%* |  1.35%* |  4.05%* |  0%
7 |  6.57%* |  0% |  1.31%* |  0%
8 |  5.12%* |  0% |  0% |  0%
9 |  5%* |  1.25%* |  -2.5% |  3.75%*
10 |  2.43%* |  1.21* |  -1.21% |  2.43%*
11 |  0% |  0% |  0% |  0%
TABLE IV: Relative improvement in fault detection effectiveness of the data diversity based approaches compared to random baseline for Linear Regression MR Set-Size | Anomaly-based | Clustering-based | Data Distribution
---|---|---|---
1 | 30.43%* | 30.43%* | 30.43%*
2 | 3.38%* | 6.77%* | 8.47%*
3 | -6.15% | 1.53%* | 9.23%*
4 | 5.79%* | 1.44%* | 5.79%*
5 | 8.33%* | 1.38%* | 13.88%*
6 | 6.75%* | 1.35%* | 10.81%*
7 | 7.89%* | 0%* | 7.89%*
8 | 6.49%* | 2.59%* | 6.49%*
9 | 5.12%* | 5.12%* | 5.12%*
10 | 2.5%* | 2.5%* | 2.5%*
11 | 0% | 0% | 0%
TABLE V: Relative improvement in fault detection effectiveness of the data diversity-based approach compared to random baseline for Naive Bayes MR Set-Size |  Rule-based |  Anomaly-based |  Clustering-based |  Data Distribution
---|---|---|---|---
1 |  38.02%* |  40.84%* |  40.84%* |  40.84*
2 |  8.69%* |  8.69%* |  8.69%* |  8.69*
3 |  2.04%* |  2.04%* |  2.04%* |  2.04%*
4 |  1.01%* |  1.01%* |  1.01%* |  1.01*
5 |  1.01%* |  1.01%* |  1.01%* |  1.01*
6 |  0% |  0% |  0% |  0%
7 |  0% |  0% |  0% |  0%
8 |  0% |  0% |  0% |  0%
9 |  0% |  0% |  0% |  0%
10 |  0% |  0% |  0% |  0%
TABLE VI: Relative improvement in fault detection effectiveness of the data diversity-based approach compared to random baseline for Convolutional Neural Network MR Set-Size | Anomaly-based | Clustering-based
---|---|---
1 | 2.38%* | 66.66%*
2 | 1.81%* | 27.27%*
3 | 1.63%* | 14.75%*
4 | 1.51%* | 6.06%*
5 | 1.44%* | 1.44%*
6 | 1.38%* | -2.77%*
7 | 1.35%* | 8.18%*
8 | 1.31%* | 5.26%*
9 | 1.29%* | 3.89%*
10 | 1.28%* | 2.56%*
11 | 0% | 0%
TABLE VII: Average relative improvement in fault detection effectiveness of data diversity-based approach over statement coverage based approach for IBk MR Set-Size |  Rule-based |  Anomaly-based |  Clustering-based |  Data Distribution
---|---|---|---|---
1 |  39.02%* |  3.20%* |  34.14%* |  14.63*
2 |  36.73%* |  26.53%* |  22.44%* |  6.12%*
3 |  9.37%* |  1.56%* |  7.81%* |  10.93%*
4 |  7.46%* |  2.98%* |  5.97%* |  1.49%*
5 |  8.69%* |  7.24%* |  10.14%* |  0%
6 |  8.33%* |  4.16%* |  6.94%* |  2.77%*
7 |  9.45%* |  2.70%* |  4.05%* |  2.70%*
8 |  7.89%* |  2.63%* |  2.63%* |  2.63%*
9 |  5.0%* |  1.25%* |  -2.50%* |  3.75%*
10 |  2.43%* |  1.21% |  -1.21% |  2.43%*
11 |  0% |  0% |  0% |  0%
TABLE VIII: Average relative improvement in fault detection effectiveness of data diversity-based approach over branch coverage based approach for IBk MR Set-Size |  Rule-based |  Anomaly-based |  Clustering-based |  Data Distribution
---|---|---|---|---
1 |  26.66%* |  26.66%* |  22.22* |  22.22%*
2 |  36.73%* |  26.53%* |  22.44%* |  6.12%*
3 |  11.11%* |  3.17%* |  9.52%* |  9.52%*
4 |  5.88%* |  1.47%* |  4.41%* |  2.49%*
5 |  5.63%* |  4.22* |  7.04%* |  2.81%*
6 |  6.84%* |  2.73%* |  5.47%* |  1.36%*
7 |  8%* |  1.33%* |  2.66%* |  1.33%
8 |  6.49%* |  1.29%* |  1.29%* |  1.29%
9 |  6.32%* |  2.53%* |  -1.26% |  5.06%*
10 |  2.43%* |  1.21%* |  -1.21%* |  2.43%*
11% |  0% |  0% |  0% |  0%
TABLE IX: Average relative improvement in fault detection effectiveness of data diversity-based approach over statement and branch coverage based approach for Naive Bayes MR Set-Size |  Rule-based |  Anomaly-based |  Clustering-based |  Data Distribution
---|---|---|---|---
1 |  40%* |  42.85%* |  42.85* |  42.85%*
2 |  9.80%* |  9.89%* |  9.89%* |  9.898%*
3 |  2.04%* |  2.04%* |  2.04%* |  2.04%*
4 |  1.01% |  1.01% |  1.01% |  1.01%
5 |  0% |  0% |  0% |  0%
6 |  0% |  0% |  0% |  0%
7 |  0% |  0% |  0% |  0%
8 |  0% |  0% |  0% |  0%
9 |  0% |  0% |  0% |  0%
10 |  0% |  0% |  0% |  0%
TABLE X: Average relative improvement in fault detection effectiveness of data diversity-based approach over statement and branch coverage based approach for Linear Regression MR Set-Size |  Anomaly-based |  Clustering-based |  Data Distribution
---|---|---|---
1 |  33.33%* |  30.43%* |  33.33%*
2 |  6.71%* |  6.77%* |  8.47%*
3 |  1.53%* |  1.53%* |  9.23%*
4 |  1.44%* |  1.44%* |  5.79%*
5 |  1.38%* |  1.38%* |  13.88%*
6 |  1.35%* |  1.35%* |  10.81%*
7 |  0% |  0% |  7.89%*
8 |  2.59%* |  2.59%* |  6.49%*
9 |  3.79%* |  5.12%* |  3.79%*
10 |  2.5%* |  2.5%* |  2.58%*
11 |  0% |  0% |  0%
TABLE XI: Average relative improvement in fault detection effectiveness of data diversity-based approach over statement and branch coverage based approach for Convolutional Neural Network MR Set-Size | Anomaly-based | Clustering-based
---|---|---
1 | 86.95%* | 204.38%*
2 | 1.81%* | 27.27%*
3 | 0% | 12.90%*
4 | 1.51%* | 6.06%*
5 | 0% | 0%
6 | 1.38%* | -2.77%
7 | 1.35%* | 8.10%*
8 | 1.31%* | 5.26%*
9 | 0% | 2.56%*
10 | 0% | 1.26%*
11 | 0% | 0%
TABLE XII: Average percentage of fault detected for subject programs Subject |  Data Distribution |  Rule Based |  Clustering Based |  Anomaly Based |  Random Based
---|---|---|---|---|---
IBk |  0.85 |  - |  0.91 |  0.72 |  0.72
Linear Regression |  0.89 |  0.92 |  0.97 |  0.95 |  0.78
CNN |  - |  - |  0.95 |  0.95 |  0.85
Naive Bayes |  0.90 |  0.93 |  0.95 |  0.92 |  0.67

####  IV-A2 Effective number of MRs used for Testing

In Figure 2, we show the effective number of MRs using the fault detection
thresholds for the subject programs. The pink vertical lines represent
effective MR set size when using the rule-based, green vertical line
represents data distribution, brown vertical line represents the clustering-
based and the orange vertical lines represent the same for the random
baseline. The respective thresholds are shown near each vertical line.

Similarly for CNN, clustering based approach provides a effective MR set size
of 2 for the 2.5% threshold and random approach provides a effective MR set
size of 4 for the 2.5% threshold. In Naive Bayes, anomaly based approach
provides a effective MR set size of 1 for the 2.5% threshold and random
approach provides a effective MR set size of 3 for the 2.5% threshold. In
Linear Regression, data diversity based approach provides an effective MR set
size of 1 for the 2.5% threshold and random approach provides a effective MR
set size of 7 for the 2.5% threshold. This practically leads to 85% reduction
in MR set size when compared to random approach. Therefore for all the subject
programs, using our proposed methods resulted in reduction in the effective MR
set size compared to using a random ordering of MRs.

####  IV-A3 The average time taken to detect a fault

Tables XIII and XIV shows the average time taken to detect a fault in Linear
Regression and IBk respectively. As shown in these results, using Data
Distribution resulted in significant reductions in the average time taken to
detect a fault by 26% when compared to random baseline. Table XV shows the
time taken to generate the metrics for the SUT. We observed that the time
taken to generate the metrics for the SUT was less than a minute. Therefore,
our proposed approach is time and cost effective when prioritizing MRs.

####  IV-A4 Average percentage of fault detected

Table XII shows the APFD for the SUT. We can observe that for all our subject
programs, our proposed metrics provided higher APFD value between 0.72-0.92
when compared to random based prioritization. As a result, our metrics
provides higher rate of fault detection when compared to random approach.

###  IV-B RQ2: Comparison of proposed MR Prioritization approach against code
coverage-based approach

In this research question, we evaluate whether the data diversity based
approach outperforms code coverage-based approach. In order to answer the
research question, we formulate the following hypotheses:

  * $H_{2}$:

For a given MR set of size $m$, the fault detection effectiveness of the MR
set produced using data diversity based approach is higher than the fault
detection effectiveness of the MR set produced by the code coverage-based
prioritization.

The null hypothesis $H_{0x}$ for each of the above defined Hypotheses $H_{x}$
is that the MR sets generated by the proposed method performs equal to or
worse than the MR sets generated by the code coverage-based prioritization in
terms of fault detection effectiveness.

Table VII and Table VIII shows the relative improvement in fault detection
percentage between the MR sets generated by data diversity based approach over
statement and branch coverage based approach for IBk. We observed that the
rule based, anomaly and clustering based metric provided improvement in fault
detection over statement coverage based approach for all the MR set sizes
except $m=11$. Therefore we can reject the null hypotheses $H_{02}$.

Table IX shows the relative improvement in fault detection percentage between
the MR sets generated by data diversity based approach and code coverage based
approach for Naive Bayes. We observed that all the metric provided improvement
in fault detection for MR set size $m=1$ to $m=4$. Therefore we cannot reject
the null hypotheses $H_{02}$ in general for Naive Bayes. Table X shows the
relative improvement in fault detection percentage between the MR sets
generated by data diversity based approach over statement and branch coverage
based approach for Linear Regression. We observed that anomaly based and data
distribution metric provided improvement in fault detection only for all MR
set sizes except $m=11$. Therefore we reject the null hypotheses $H_{02}$ in
general for Linear Regression. Table XI shows the relative improvement in
fault detection percentage between the MR sets generated by data diversity
based approach over statement and branch coverage based approach for CNN. We
observed that the clustering-based approach showed improvement in fault
detection for all the MR set sizes and anomaly based approach showed
improvement in fault detection only for MR set size $m=1$ to $m=5$. Therefore,
we reject the null hypothesis $H_{02}$ in general for clustering-based metric.
Overall, our proposed metrics outperforms the statement and branch coverage
based approach for the subject programs.

TABLE XIII: Comparison of average time taken to detect a fault using Data Diversity method and Random baseline for Linear Regression and IBk Subject |  Tsp |  Data Distribution |  Rule Based |  Random Baseline |  % Time Reduction
---|---|---|---|---|---
IBk |  Test1 |  - |  8910.80s |  12562.49s |  29.06%
Linear Regression |  Prediction |  17057.43s |  - |  21609.766s |  26.06%
TABLE XIV: Comparison of average time taken to detect a fault using Data Diversity method and Random baseline for IBk $\mathbf{{T_{sv}}}$ |  $\mathbf{{F_{v}}}$ |  Rule Based |  Random baseline |  % Time reduction
---|---|---|---|---
Testset1 |  PIT |  8910.80s |  12562.49s |  29.06%
Testset2 |  Major |  6148.29s |  7508.29s |  22.23%
TABLE XV: Time taken to generate the metrics for SUT SUT |  Rule-based |  Anomaly-based |  Clustering-based |  Data Distribution
---|---|---|---|---
IBk |  9 sec |  7 sec |  17 sec |  6 sec
Linear Regression |  - |  4 sec |  20 sec |  3 sec
Naive Bayes |  8 sec |  6 sec |  18 sec |  5 sec
CNN |  N/A |  200 sec |  300 sec |  N/A

##  V Related Work

Prior work [1] proposed fault-based and code coverage-based approach to
prioritize MRs in the context of regression testing. The fault-based approach
applies a greedy algorithm on the faults detected by the MRs on the previous
version of the SUT to prioritize MRs for testing the next version of the SUT.
Code coverage-based approach uses statements and branches covered by the MRs
to prioritize MRs. The proposed approach was applied to three diverse
applications and the result indicates that the code coverage based approach
did not work effectively for machine learning-based programs.

Mayer and Guderlei [28] conducted an empirical study to evaluate the
usefulness of MRs. Based on the experiment result, the authors devise several
criteria to asses the quality of the MRs based on the potential usefulness.
These criteria were defined based on the different execution paths undertaken
by the source and follow-up test cases. While this work provides some
directions into selecting useful MRs, they do not directly address MR
prioritization. In particular, we focus on developing automated methods for MR
prioritization. Hui et al. [29] provided qualitative guidelines for MR
prioritization. The guidelines could be used to select effective MRs and
improve the performance of MT. The author also provides three measurable
metrics for MR selection: in-degree of MR, algebra complexity of MR, and
distance between test inputs of MR. The empirical result shows inconsistency
between the correlation value and the metrics, indicating that other factors
could influence the effectiveness of MR. Mani et al. [30] proposed metrics to
identify the quality of dataset based on centroid positioning, pairwise
boundary conditioning and equivalence partitioning. This work does not provide
insight on the diversity in the dataset and does not address MR
prioritization. Sun et al. [31] proposed a technique to generate good source
test cases in metamorphic testing based on constraint solvers and techniques
of symbolic execution. Also, a prioritization of source test cases was
conducted based on path distances among test cases. This work does not focus
on prioritization of metamorphic relations. Also, the work does not conduct
prioritization of metamorphic relations for machine learning programs.

##  VI Threats to Validity

External Validity:  Our proposed metric is applied to supervised ML
classifiers. Metrics such as rule-based is not suitable for unsupervised
classifiers and data-distribution metric can be applied only to data set with
numerical values. We plan to design more metrics to cover unsupervised
classifiers and non-numerical data. In this work, all the experiment subjects
are implementation of machine learning algorithm used in the Weka ML library.
Although they are popular algorithms and widely used, our findings may not be
generalizable to commercial machine learning projects such as autonomous
vehicles. To mitigate this threat, we plan to explore the effectiveness of our
proposed approach on AI systems in the health care domain and Apollo
autonomous driving platform projects in the future.

Internal Validity:  We used dataset having an average of 500000 instances and
90 attributes on our subject programs. However, the dataset size could be
small for generating the proposed metrics and prioritizing the MRs. We used
less than 100 mutants for validating our prioritized MR ordering for CNN
program. However, it is possible that the number of mutants used is low.

##  VII Conclusion

Metamorphic relations in metamorphic testing have multiple source and follow-
up test cases. As a result, execution and training of machine learning model
with large source and follow-up data set can exponentially increase the
execution time and cost of metamorphic testing for ML applications. To
overcome this problem, we developed four metrics based on diversity between
the source and follow-up data set for prioritization of MRs.

We evaluated our proposed metrics on four open-source machine learning
programs. The experimental results show that utilizing MR prioritization on ML
applications can increase the fault detection effectiveness of a given MR set
up to 33% for linear regression, 86% for CNN and up to 40% for Naive Bayes. In
particular, using our proposed MR prioritization approaches reduced the
average time taken to detect a fault compared to most common approach of
randomly executing the MRs by 30% for KNN classifier implementation. Further,
our results show that data diversity based metrics outperforms code coverage
based approach for prioritizing MRs by 40%. Our proposed metrics also provided
higher rate of fault detection when compared to random based prioritization.
Our finding also indicates that the number of MRs needed to execute and test
ML programs reduced by 66% which translates to saving cost and time for
software testing practitioners.

## References

  * [1] Anonymous, “Under submission. details omitted for double-blind reviewing.”
  * [2] S. Liao, “Chinese facial recognition system mistakes a face on a bus for a jaywalker.”
  * [3] C. Ziegler, “A google self-driving car caused a crash for the first time,” 2016\.
  * [4] J. Glenday, “Ai-powered genderify platform shut down after bias-based backlash.”
  * [5] A. Ohnsman, “Lidar maker velodyne ‘baffled’by self-driving uber’s failure to avoid pedestrian,” 2018.
  * [6] E. T. Barr, M. Harman, P. McMinn, M. Shahbaz, and S. Yoo, “The oracle problem in software testing: A survey,” _IEEE transactions on software engineering_ , vol. 41, no. 5, pp. 507–525, 2015.
  * [7] E. J. Weyuker, “On testing non-testable programs,” _The Computer Journal_ , vol. 25, no. 4, pp. 465–470, 1982.
  * [8] J. M. Zhang, M. Harman, L. Ma, and Y. Liu, “Machine learning testing: Survey, landscapes and horizons,” _IEEE Transactions on Software Engineering_ , 2020\.
  * [9] T. Y. Chen, F.-C. Kuo, H. Liu, P.-L. Poon, D. Towey, T. Tse, and Z. Q. Zhou, “Metamorphic testing: A review of challenges and opportunities,” _ACM Computing Surveys (CSUR)_ , vol. 51, no. 1, p. 4, 2018.
  * [10] S. Segura and Z. Q. Zhou, “Metamorphic testing 20 years later: A hands-on introduction,” in _Proceedings of the 40th International Conference on Software Engineering: Companion Proceeedings_. ACM, 2018, pp. 538–539.
  * [11] C. Murphy, G. E. Kaiser, and L. Hu, “Properties of machine learning applications for use in metamorphic testing,” 2008.
  * [12] X. Xie, J. W. Ho, C. Murphy, G. Kaiser, B. Xu, and T. Y. Chen, “Testing and validating machine learning classifiers by metamorphic testing,” _Journal of Systems and Software_ , vol. 84, no. 4, pp. 544–558, 2011.
  * [13] X. Xie, J. Ho, C. Murphy, G. Kaiser, B. Xu, and T. Y. Chen, “Application of metamorphic testing to supervised classifiers,” in _2009 Ninth International Conference on Quality Software_. IEEE, 2009, pp. 135–144.
  * [14] J. E. Gewehr, M. Szugat, and R. Zimmer, “Bioweka—extending the weka framework for bioinformatics,” _Bioinformatics_ , vol. 23, no. 5, pp. 651–653, 2007.
  * [15] M. Zhang, Y. Zhang, L. Zhang, C. Liu, and S. Khurshid, “Deeproad: Gan-based metamorphic testing and input validation framework for autonomous driving systems,” in _2018 33rd IEEE/ACM International Conference on Automated Software Engineering (ASE)_. IEEE, 2018, pp. 132–142.
  * [16] Anonymous, “Under submission. details omitted for double-blind reviewing.”
  * [17] G. Gay, M. Staats, M. Whalen, and M. P. Heimdahl, “Automated oracle data selection support,” _IEEE Transactions on Software Engineering_ , vol. 41, no. 11, pp. 1119–1137, 2015.
  * [18] S. Elbaum, A. G. Malishevsky, and G. Rothermel, “Test case prioritization: A family of empirical studies,” _IEEE transactions on software engineering_ , vol. 28, no. 2, pp. 159–182, 2002.
  * [19] A. G. Malishevsky, J. R. Ruthruff, G. Rothermel, and S. Elbaum, “Cost-cognizant test case prioritization,” Citeseer, Tech. Rep., 2006.
  * [20] S. Elbaum, G. Rothermel, S. Kanduri, and A. G. Malishevsky, “Selecting a cost-effective test case prioritization technique,” _Software Quality Journal_ , vol. 12, no. 3, pp. 185–210, 2004.
  * [21] D. W. Aha, D. Kibler, and M. K. Albert, “Instance-based learning algorithms,” _Machine Learning_ , vol. 6, no. 1, pp. 37–66, Jan 1991. [Online]. Available: https://doi.org/10.1007/BF00153759
  * [22] Q.-H. Luu, M. F. Lau, S. P. Ng, and T. Y. Chen, “Testing multiple linear regression systems with metamorphic testing,” _Journal of Systems and Software_ , vol. 182, p. 111062, 2021.
  * [23] J. Ding, X. Kang, and X.-H. Hu, “Validating a deep learning framework by metamorphic testing,” in _2017 IEEE/ACM 2nd International Workshop on Metamorphic Testing (MET)_. IEEE, 2017, pp. 28–34.
  * [24] Y.-S. Ma, J. Offutt, and Y.-R. Kwon, “Mujava: a mutation system for java,” in _Proceedings of the 28th international conference on Software engineering_. ACM, 2006, pp. 827–830.
  * [25] H. Coles, T. Laurent, C. Henard, M. Papadakis, and A. Ventresque, “Pit: a practical mutation testing tool for java,” in _Proceedings of the 25th International Symposium on Software Testing and Analysis_. ACM, 2016, pp. 449–452.
  * [26] R. Just, “The major mutation framework: Efficient and scalable mutation analysis for java,” in _Proceedings of the 2014 International Symposium on Software Testing and Analysis_. ACM, 2014, pp. 433–436.
  * [27] P. Good, _Permutation tests: a practical guide to resampling methods for testing hypotheses_. Springer Science & Business Media, 2013.
  * [28] J. Mayer and R. Guderlei, “An empirical study on the selection of good metamorphic relations,” in _30th Annual International Computer Software and Applications Conference (COMPSAC’06)_ , vol. 1. IEEE, 2006, pp. 475–484.
  * [29] Z.-W. Hui, S. Huang, H. Li, J.-H. Liu, and L.-P. Rao, “Measurable metrics for qualitative guidelines of metamorphic relation,” in _2015 IEEE 39th Annual Computer Software and Applications Conference_ , vol. 3. IEEE, 2015, pp. 417–422.
  * [30] S. Mani, A. Sankaran, S. Tamilselvam, and A. Sethi, “Coverage testing of deep learning models using dataset characterization,” _arXiv preprint arXiv:1911.07309_ , 2019.
  * [31] “Path-directed source test case generation and prioritization in metamorphic testing,” _Journal of Systems and Software_ , vol. 183, p. 111091, 2022. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0164121221001886