# BlockEcho: Retaining Long-Range Dependencies for Imputing Block-Wise Missing
Data

Qiao Han  Mingqian Li  Yao Yang&Yiteng Zhai∗
Zhejiang Lab
{hanq, mingqian.li, yangyao<EMAIL_ADDRESS>
###### Abstract

Block-wise missing data poses significant challenges in real-world data
imputation tasks. Compared to scattered missing data, block-wise gaps
exacerbate adverse effects on subsequent analytic and machine learning tasks,
as the lack of local neighboring elements significantly reduces the
interpolation capability and predictive power. However, this issue has not
received adequate attention. Most SOTA matrix completion methods appeared less
effective, primarily due to overreliance on neighboring elements for
predictions. We systematically analyze the issue and propose a novel matrix
completion method “BlockEcho” for a more comprehensive solution. This method
creatively integrates Matrix Factorization (MF) within Generative Adversarial
Networks (GAN) to explicitly retain long-distance inter-element relationships
in the original matrix. Besides, we incorporate an additional discriminator
for GAN, comparing the generator’s intermediate progress with pre-trained MF
results to constrain high-order feature distributions. Subsequently, we
evaluate BlockEcho on public datasets across three domains. Results
demonstrate superior performance over both traditional and SOTA methods when
imputing block-wise missing data, especially at higher missing rates. The
advantage also holds for scattered missing data at high missing rates. We also
contribute on the analyses in providing theoretical justification on the
optimality and convergence of fusing MF and GAN for missing block data.

##  1 Introduction

The issue of missing data is prevalent in real-world datasets, stemming from
factors such as users’ reluctance to share data, variations in feature
structures, privacy concerns, and unintended data corruption over time
Scheffer (2002). In scenarios involving large-scale and highly correlated
missing data, matrix completion becomes essential to impute the unknown
entries accurately and efficiently. Accurately and efficiently estimating the
unknown elements in a matrix is pivotal for effective dataset analysis and
supports subsequent machine learning and operations research tasks. Matrix
completion finds wide application in recommender systems, trajectory recovery,
phase retrieval, computer vision, genotype imputation and other fields of
research Santos et al. (2019). Missing data isn’t limited to randomly
scattered elements; it can also occur in blocks, referred to as block-wise
missing data. For example, discrepancies in feature sets provided by different
agents can result in non-overlapping blocks of NA values during data
aggregation. Moreover, certain regions may experience skipped clinical trials
due to technical and financial constraints, while continuous data corruption
can arise from malfunctioning monitoring devices over time.

A diverse array of methods have been developed over the years to address the
matrix completion problem, spanning traditional techniques like linear
interpolation, K-nearest neighbors imputation, Multiple Imputation by Chained
Equations (MICE), and ensemble methods such as MissForest Williams (2015). In
the new data regime, Generative Adversarial Networks (GANs) have also
demonstrated immense potential in accurately estimating high-dimensional data
distributions Goodfellow et al. (2020). Consequently, recent works have
explored GAN architectures for data recovery in the missing data context.
However, most prevailing completion techniques display suboptimal
effectiveness on block-missing matrices, especially with increasing missing
proportions. They predominantly leverage neighboring available elements to
predict unknown entries, making them vulnerable to systematic feature gaps. An
exception is Matrix Factorization (MF) which implicitly retains inter-sample
dependencies across longer ranges Hastie et al. (2015). But its linearity
assumptions limit capturing complex nonlinear relationships. To unite their
complementary strengths and offset limitations, we propose BlockEcho \- an
integrated approach tailored for block-missing data that capitalizes on the
strengths of both GAN and MF. It explicitly retains long-range inter-element
associations via MF while modeling intrinsic nonlinear patterns through GANs.
The key contributions are:

  1. 1.

We have conducted an in-depth analysis and formally defined the concept of
“Block-wise” missing data. Subsequently, we have innovatively introduced a
solution to address this challenge, named BlockEcho.

  2. 2.

We extensively benchmark BlockEcho against SOTA baselines on diverse public
datasets for traffic, epidemiology, and recommendations. Results demonstrate
significant gains in imputation accuracy under both block-missing and high-
rate scattered missing regimes. Downstream forecasting tasks further showcase
advantages.

  3. 3.

We provide theoretical analysis into global optimality and convergence for the
proposed integrated objective. We particularly justify the synergistic effects
of MF encoding long-range dependencies and GANs locally adapting complex
distributions. Rigorous proofs supplement empirical observations.

The rest of this paper is structured as follows. In Section 2 we provide a
concise literature review in missing data imputation studies. Section 3
defines the problem. We elaborate our model BlockEcho in Section 4. Section 5
provides a theoretical discussion. Section 6 details our experiments and
results. We conclude the paper in Section 7.

Figure 1: Problem definition, obejective and evaluation metric

##  2 Related Work

###  2.1 Data imputation approaches

Broadly, prevailing techniques for missing data imputation fall under two
categories – discriminative or generative. Discriminative methods directly
estimate the conditional distribution of missing values given observed
entries. This includes MICE Buuren and Oudshoorn (2000), MissForest Stekhoven
and Bühlmann (2011) and Matrix Factorization Hastie et al. (2015) Yu et al.
(2016). Generative approaches instead model the complete joint data
distribution and draw imputes from conditionals on observed parts.
Traditionally, expectation maximization algorithms iteratively refining
estimates were common Gold and Bentler (2000). Recently, advances in deep
generative modeling have spawned powerful imputation paradigms. Notably, J.
Yoon et al. Yoon et al. (2018) pioneered Generative Adversarial Networks
(GANs) for missing data imputation through the GAIN framework in 2018. Instead
of assuming explicit parametric forms, GANs can flexibly approximate arbitrary
data distributions. Consequently, GAN-based solutions have risen as SOTA,
significantly improving over pre-GAN techniques.

###  2.2 GAN-based methods

Many works have since built upon GAIN, either enhancing its convergence like
SGAIN (slim Wasserstein GAIN) Neves et al. (2021), integrating auxiliary
models as in GAMIN (Generative adversarial multiple imputation network) Yoon
and Sull (2020), or exploiting complementary data modalities like images via
CollaGAN (Collaborative GAN for missing image data imputation) Lee et al.
(2019). Some explore sophisticated relationships within data beyond
distributions, including: 1) Dual generative modeling of values and high-
dimensional feature abstractions in MisGAN Li et al. (2019), 2) time-series
self-attention mechanisms in Zhang et al. (2021) to capture sensor inter-
dependencies, 3) Pretraining classifier-guided generative models as PC-GAIN
Wang et al. (2021), 4) Distribution-centered losses like STGAN Yuan et al.
(2022), and 5) Multimodal spatio-temporal modeling via GAN+RNN+GCN Kong et al.
(2023). ANODE-GAN Chang et al. (2023) also augments GANs with auxiliary
variational autoencoders.

Nonetheless, prevailing GAN imputation approaches overlook systematic data
deficiencies like block-missing data. Mostly relying on local neighborhoods,
they remain less effective as gaps exacerbate, especially at higher missing
rates. Our proposed BlockEcho framework aims to address this limitation by
uniquely blending GANs with matrix factorization.

##  3 Problem Definition

Consider an original data matrix $X\in\mathbb{R}^{m*n}$. $M\in\\{0,1\\}^{m*n}$
is set as a binary mask matrix codifying missingness, where $M_{ij}=0$ denotes
element $X_{ij}$ as unobserved. Imposing $M$ on $X$ gives the incomplete
matrix:

$\tilde{X}_{ij}=\begin{cases}X_{ij}&ifM_{ij}=1\\\ NaN&ifM_{ij}=0\\\
\end{cases}$ (1)

We refer to $X$ as the original matrix, $M$ as the mask matrix, and
$\tilde{X}$ as the missing data matrix. The data imputation problem aims to
obtain an estimation $\hat{X}$ of $X$ given $\tilde{X}$ and $M$, with the
objective of minimizing the discrepancy between the imputed data and the
original data. Our final output is the imputed matrix $\bar{X}$ defined as

$\bar{X}:=\tilde{X}\odot M+\hat{X}\odot(1-M)$ (2)

where $\odot$ denotes element-wise multiplication.

For any given coordinates $(i_{l},j_{l})$, $\exists\ i_{u}\geq i_{l}+3$ and
$j_{u}\geq j_{l}+3$, $\forall\ i_{l}\leq i\leq i_{u}$ and $j_{l}\leq j\leq
j_{u}$, $M_{ij}=0$, then we define this area as exhibiting a local block-wise
missing pattern. Following the definition, the missing data in datasets can
manifest as uni-block-wise, multiple-separate-block-wise, scattered, or as a
mixed combination thereof.

We conducted a straightforward experiment to understand the harm of missing
data. We employed a random forest algorithm to forecast the next timestamp’s
traffic flow for dataset PE-BAY (outlined in Section 7), without any feature
engineering but normalization. Each category of missing data exhibited a 60%
missing rate. Our experimental findings, as presented in Table 1, highlighted
the notably adverse impact of a substantial, uniformly absent data block on
the accuracy of predictions. Henceforth, unless explicitly specified, “block-
wise missing” denotes the uni-block.

Ori-Data | Scattered | Tri-Blocks | Uni-Block
---|---|---|---
0.0395 | 0.0961 | 0.1132 | 0.1249
Table 1: WMAPE of prediction task in different missing type  Figure 2: The
architecture of BlockEcho

##  4 The Architecture of “BlockEcho”

In this section, we will elaborate on our integrated model, elucidating how
the MF and GAN components synergize to enhance overall performance. The
overarching model architecture of BlockEcho is depicted in Figure 2.

###  4.1 Matrix Factorization (MF) and Pre-training

Within the MF component, the missing data matrix is decomposed into two
smaller matrices, denoted as the imputed row embedding matrix $U$ and the
imputed column embedding matrix $V$. Their dimensions are $m*h$ (adjustable to
mini_batch_size$*h$ during training) and $h*n$, where $h$ is a hyperparameter.
In essence, $h$ serves as an indicator of the original matrix’s level of
freedom, with smaller values reflecting a higher degree of interdependency
among the data in the original matrix.

As illustrated in Figure 2, there are two MF components in BlockEcho. One is
deployed during the pre-training stage, breaking down the matrix into $U_{p}$
and $V_{p}$ to extract long-range relationship features, which are then used
to supervise model training in a transformed space through a discriminator
$D_{I}$. The other is directly involved in neural network training and result
generation: the imputed row embedding matrix $U$ is produced by the generator,
while the imputed column embedding matrix $V$ is directly calculated through
backpropagation. Subsequently, these matrices undergo a matrix completion
layer (MCL) to match the known parts of the missing data matrix. Rather than
employing direct multiplication, MCL incorporates a fully connected neural
network layer after multiplication. . This technique is employed to introduce
nonlinear relationships while upholding the inherent advantage of low
constraints on the degree of freedom in the matrix factorization method. Both
MF components are oriented towards minimizing the Kullback-Leibler divergence
of known elements, thereby establishing the MF objective as

$\begin{split}\min\ \ &\sum_{i,j}((X\odot M)_{ij}log(\frac{(X\odot
M)_{ij}}{(MCL(UV)\odot M)_{ij}})\\\ &-(X\odot M)_{ij}+(MCL(UV)\odot
M)_{ij})\end{split}$ (3)

###  4.2 Generative Adversarial Nets (GANs)

The GAN framework consists of a generator $G$ and two discriminators
$D_{I},D_{II}$ as denoted in Figure 2. The generator $G$ takes $\tilde{X}$
with 0 imputed, $M$, and a noise matrix $Z$ with the same size as $U$ as the
inputs, and $U$ as the output. Next, MCL takes $U$ as input and outputs the
estimation $\hat{X}$. With $\hat{X}$, $\tilde{X}$ and $M$, we build the
imputed matrix $\bar{X}$ following Equation (2). Given a vector of length h,
the discriminator $D_{I}$ will determine whether it originates from pre-
trained matrix factorization or from the generator. To do this, we randomly
generate $Y\in\\{0,1\\}^{m*1}$ as the index vector and build a mixed-row
matrix $U_{D}$ combining $U_{p}$ and $U$:

$U_{D}=(YI)\odot U_{p}+(1-YI)\odot U$ (4)

where $I\in\\{1\\}^{1*h}$. The discriminator $D_{I}$ will output a result of
the same size as $Y$ and the objective is

$\mathop{\min}\limits_{G}\mathop{\max}\limits_{D_{I}}\
Y^{T}logD_{I}(U_{D})+(1-Y)^{T}log(1-D_{I}(U_{D}))$ (5)

The discriminator $D_{II}$ is designed to discern whether the elements in
$\bar{X}$ is from $\tilde{X}$ (real) or $\hat{X}$ (fake). This differs from a
standard GAN model which aims to differentiate inputs ranging from entirely
true to entirely false. In addition to $\bar{X}$, $D_{II}$ incorporates the
hint matrix $H$ as input. Hint matrix, proposed in GAIN method Yoon et al.
(2018), modifies the elements of $M$ to 1/2 at a specified rate. This
modification facilitates $D_{II}$ to converge more swiftly and accurately by
furnishing “enough” information regarding $M$. The objective of this component
is

$\begin{split}\mathop{\min}\limits_{G}\mathop{\max}\limits_{D_{II}}&\
\sum_{i,j}(M\odot logD_{II}(\bar{X},H)\\\ &+(1-M)\odot
log(1-D_{II}(\bar{X},H)))_{ij}\end{split}$ (6)

In particular, the MF component serves as a structured constraint and guidance
for the GAN component (analogous to the relationship between CNN and MLP),
enabling the entire network to explicitly preserve long-range relationships
(similar to CNN’s explicit preservation of proximity relationships). This
fundamental characteristic underpins the superior performance of our model
compared to baselines on datasets with missing blocks or high missing rates.

###  4.3 Objective

According to (3), (5) and (6), the objective of the whole multi-loss model is

$\begin{split}\mathop{\min}\limits_{G}\mathop{\max}\limits_{D}&\ \
(1-\alpha)[(1-Y)^{T}log(1-D_{I}(U_{D}))\\\
&+Y^{T}logD_{I}(U_{D})+\sum_{i,j}(M\odot logD_{II}(\bar{X},H)\\\ &+(1-M)\odot
log(1-D_{II}(\bar{X},H)))_{ij}]\\\ &+\alpha[\sum_{i,j}((X\odot
M)_{ij}log(\frac{(X\odot M)_{ij}}{(G(\tilde{X},M,Z)V\odot M)_{ij}})\\\
&-(X\odot M)_{ij}+(G(\tilde{X},M,Z)V\odot M)_{ij})]\\\ with\ \
&U_{D}=(YI)\odot U_{p}+(1-YI)\odot G(\tilde{X},M,Z),\\\
&\bar{X}=\tilde{X}\odot M+MCL(G(\tilde{X},M,Z)V)\odot(1-M)\end{split}$

(7)

Similar to a standard GAN, we alternate between training the generator and
discriminators separately. The MF loss (3) is only connected to the generator
and will be utilized for multi-loss training alongside the generator’s
objective.

##  5 Theoretical Discussion

###  5.1 Global Optimality

Without loss of generality, we assume the elements of $X$ as
$X_{ij}=f(i,j)+s_{ij}$. Here, $f$ represents the bias effect of coordinates on
X. The greater the internal correlation of $X$’s elements, the stronger
regularity $f$ shows. $\\{s_{ij}:i\in[0,m-1]\,j\in{0,n-1}\\}$ are a series of
random variables which follow independent and identical distribution $p_{d}$.
We can define $\hat{X}$ in same way, i.e.,
$\hat{X}_{ij}=\hat{f}(i,j)+\hat{s}_{ij}$, where $\hat{s}$’s probability
distribution is $p_{g}$. In real-world data, whether an element is masked or
not may be affected by the value of the element itself. It also has different
effects on the corresponding elements of $\hat{X}$ in generative model. So we
need to discuss the probability distribution under different mask conditions.
We assume that the distributions of $s$ corresponding to the unmasked elements
in $X$ and $\hat{X}$ are $p_{d|m=1}$ and $p_{g|m=1}$, and to other masked
elements are $p_{d|m=0}$ and $p_{g|m=0}$. The objective of this problem can be
written as

$\hat{f}\to f\ \ and\ \ p_{g|m=0}\to p_{d|m=0}$ (8)

Throughout this section, we assume that the data are at MCAR (missing
completely at random)Heitjan and Basu (1996), i.e. $M$ is independent of $X$.
Then we can obtain $p_{d|m=0}=p_{d|m=1}=p_{d}$, and (8) can be transformed
into

$\hat{f}\to f\ \ and\ \ p_{g|m=0}\to p_{d|m=1},$ (9)

or

$\hat{f}\to f\ \ and\ \ p_{g|m=1}\to p_{d|m=1}\ \ and\ \ p_{g|m=0}\to
p_{g|m=1}.$ (10)

Obviously (8), (9) and (10) are equivalent.

Below we will prove that the theoretical global optimal of (6) is the same as
that in (9). The loss function for standard GAN is

$E_{x\sim p_{d}}[logD(x)]+E_{z\sim p_{g}}[log(1-D(G(z)))]$ (11)

For this question, it can be re-written as:

$\begin{split}&E_{\bar{S}\sim p_{d|m=1}}[logD(\bar{X},H)]+E_{\bar{S}\sim
p_{g|m=0}}[log(1-D(\bar{X},H))]\\\ &=\int_{x}p_{d|m=1}(x-f)logD(x,H)\\\
&+p_{g|m=0}(x-\hat{f})log(1-D(x,H))\mathrm{d}x\end{split}$ (12)

For any $a,b>0$, the function $alogx+blog(1-x)$ reaches its minimum if and
only if $x=\frac{a}{a+b}$. So if G is fixed, the optimal discriminator D is
$\frac{p_{d|m=1}(x-f)}{p_{d|m=1}(x-f)+p_{g|m=0}(x-\hat{f})}$. Substituting the
optimal discriminator into (12), we can obtain

$\begin{aligned} &E_{\bar{S}\sim p_{d|m=1}}[logD(\bar{X},H)]+E_{\bar{S}\sim
p_{g|m=0}}[log(1-D(\bar{X},H))]\\\
&=KL(p_{d|m=1}(x-f)\|\frac{1}{2}(p_{d|m=1}(x-f)+p_{g|m=0}(x-\hat{f})))\\\
&+KL(p_{g|m=0}(x-\hat{f})\|\frac{1}{2}(p_{d|m=1}(x-f)+p_{g|m=0}(x-\hat{f})))\\\
&-log(4)\end{aligned}$

(13)

We recognize that (12) is essentially optimizing Kullback-Leibler divergence.
It will achieve its global minimum $-log(4)$ if and only if $\hat{f}=f\ \&\
p_{g|m=0}=p_{d|m=1}$, which is the same as (9).

Next, we will calculate the theoretical global optimal of (3). To do this, we
initially introduce the entropy of the variable, $Ent(A)$, a measure of the
disorder of $A$. According to the definitions of entropy and mutual
information, we can find that

$\begin{split}Ent(\hat{X})=&Ent(\hat{X}\odot M+\hat{X}\odot(1-M))\\\
=&Ent(\hat{X}\odot M)+Ent(\hat{X}\odot(1-M))\\\ &-I(\hat{X}\odot
M,\hat{X}\odot(1-M))\end{split}$ (14)

where $I$ denotes mutual information. Noted that $I(\hat{X}\odot
M,\hat{X}\odot(1-M))\leq Ent(\hat{X}\odot(1-M))$ and the equality holds when
$p_{g|m=0}=p_{g|m=1}$, which means that if $p_{g|m=1}$ is fixed, the global
optimal of $min(Ent(\hat{X}))$ is $p_{g|m=0}=p_{g|m=1}$. Therefore, (10) is
equivalent to a multi-objective optimization problem

$min\ [D(X\odot M||\hat{X}\odot M),\ Ent(\hat{X})]$ (15)

Where $D(\hat{X}\odot M||X\odot M)$ is the primary objective, which can be
Euclidean distance or Kullback-Leibler divergence. Its theoretical optimal
solution is $\hat{f}=f\ \&\ p_{g|m=1}=p_{d|m=1}$ under the assumption that
$\hat{X}$ could be fitted by any functions (actually not). $Ent(\hat{X})$ is
the second objective. We assume the global optimality of (15) is
$p_{g|m=0}=p_{g|m=1}=p_{d|m=1}=p^{*}$ , in which case
$Ent(\hat{X})=Ent^{*}(\hat{X})$. Thus, we have shown that problem (10) is
equivalent to the multi-objective optimization problem (15).

Actually, it is very difficult to quantify $Ent(\hat{X})$ in (15) in practice.
Matrix factorization is essentially a constraint on entropy in the form of the
network structure. From this premise, we target at solving (15) approximately,
and the optimization problem for the matrix factorization can be written as:

$\begin{split}min\ \ &D(X\odot M||\hat{X}\odot M)\\\ s.t.\ \
&Ent(\hat{X})<ent\end{split}$ (16)

Where $ent$ is a hyperparameter decided by the network structure (mainly $h$
mentioned in section 4.1). We proceed with the following two cases: 1\.
$ent>Ent^{*}(\hat{X})$. In this case, the optimal solution of $\hat{f}=f\ \&\
p_{g|m=1}\to p_{d|m=1}$ can be reached theoretically, but it lacks a
sufficient solution for $p_{g|m=0}\to p_{g|m=1}$. It corresponds to
overfitting in general machine learning models. 2. $ent<Ent^{*}(\hat{X})$. In
this case, $\hat{x}$ cannot be fully represented by the model due to its low
degree of freedom limited by the tighter constraint on entropy, which hinders
$\hat{x}$ from reaching the optimal of both $\hat{f}=f\ \&\ p_{g|m=1}\to
p_{d|m=1}$ and $p_{g|m=0}\to p_{g|m=1}$. It corresponds to underfitting in
general machine learning models. In either case, the theoretical global
optimal of MF, i.e., problem (16), is the approximate solution of (15).

| Method | ME-LA | PE-BAY | Cov-ca | Cov-de | Movie
---|---|---|---|---|---|---
| MissForest | 0.2166 | 0.1037 | 0.1106 | 0.1589 | –
Block-wise | Matrix Factorization | 0.1637 | 0.0814 | 0.0645 | 0.0802 | –
Missing | GAIN | 0.2020 | 0.1135 | 0.0642 | 0.0811 | –
Data | PC-GAIN | 0.2169 | 0.1340 | 0.1055 | 0.1166 | –
| STGAN | 0.2057 | 0.1229 | 0.0757 | 0.0968 | –
| BlockEcho(ours) | 0.1561 | 0.0743 | 0.0495 | 0.0784 | –
| MissForest | 0.1328 | 0.0612 | 0.0493 | 0.0584 | 0.1998
Scattered | Matrix Factorization | 0.1633 | 0.0816 | 0.0593 | 0.0737 | 0.1821
Missing | GAIN | 0.1583 | 0.0664 | 0.0471 | 0.0616 | 0.2248
Data | PC-GAIN | 0.1327 | 0.0617 | 0.0472 | 0.0673 | 0.1937
| STGAN | 0.1444 | 0.0679 | 0.0351 | 0.0470 | 0.1983
| BlockEcho(ours) | 0.1322 | 0.0608 | 0.0321 | 0.0488 | 0.1797
Table 2: RMSE performance comparison of different data imputation methods at a
fixed missing rate of 60%. Noted that in Movie dataset, 80% of the data are
inherently missing, so we choose to mask the remaining visible elements at a
missing rate of 60%.

(a) Block-wise missing data

(b) Scattered missing data

Figure 3: RMSE performance trend with increasing missing rate

###  5.2 Convergence

In the previous subsection, we demonstrated that loss (3) has a similar but
relatively inferior theoretical optimal solution compared to loss (6), even
when we introduce non-linear transformations to MF through neural networks.
However, it is anticipated intuitively (also be verified by ablation
experiments in subsection 6.3), that the long-range dependencies introduced by
MF will exert a significant influence on the outcomes of block-missing data
imputation. This is due to the inherent challenge faced by models with a
higher degree of freedom in effectively converging to fit long-range
dependencies, particularly in the case of GAN models. The convergence capacity
of GANs has posed a substantial challenge since its inception. In 2017, L
Mescheder et al. Mescheder et al. (2017) postulated that poor convergence is
attributed to the eigenvalues of the Jacobian matrix of the gradient vector
field having a zero real part and an excessively large imaginary part.
Consequently, strategies such as gradient penalty Gao et al. (2020) and
spectral normalization Miyato et al. (2018) have been adopted to alleviate
this limitation. In tackling the specific matrix completion problem at hand,
we leverage (3), the MF loss function, and (5), which enables the generator to
more accurately emulate MF, to effectively steer the convergence of (6).

At last of the section, we will prove the convergence of loss (3), indicating
that $D(\hat{X}\odot M||X\odot M)$ is nonincreasing under certain update
rules.

Let $F$ denote this expression. Let $u$ denote the set of parameters of the
model, and $u^{t}$ the set of parameters upon update at step $t$. We introduce
an auxiliary function $G(u,u^{t})$ for $F(u)$ satisfying $G(u,u)=F(u)$ and
$G(u,u^{t})\geq F(u)$.Consequently, we can obtain $F(u^{t+1})\leq F(u^{t})$
when $G(u^{t+1},u^{t})\leq G(u^{t},u^{t})$, then $F$ is nonincreasing under
the update $u^{t+1}=argmin_{u}G(u,u^{t})$.

In our paper, we choose Kullback-Leibler divergence as the loss (3) because
the datasets are closer to Poisson distributions. In 2001, Daniel et al.Lee
and Seung (2000) proposed the general auxiliary function for divergence. For
this problem, we can make $F(u)$ equal to loss (3) as
$\hat{X}_{i,j}:=\sum_{a}u_{ia}v_{aj}$ and
$\hat{X}_{i,j}^{t}:=\sum_{a}u_{ia}^{t}v_{aj}$, and define $G(u,u^{t})$ as

$\begin{split}&G(u,u^{t}):=\sum_{i,j}(((X\odot M)_{ij}log(X\odot M)_{ij}\\\
&-(X\odot M)_{ij}+(\hat{X}\odot M)_{ij})\\\ &-\sum_{a}(X\odot
M)_{ij}\frac{u_{ia}^{t}v_{aj}}{\hat{X}_{i,j}^{t}}(logu_{ia}v_{aj}-log\frac{u_{ia}^{t}v_{aj}}{(\hat{X}\odot
M)_{ij}^{t}}))\end{split}$ (17)

It is straightforward to verify that $G$ satisfies $G(u,u)=F(u)$ and
$G(u,u^{t})>=F(u)$ , and that $F(u^{t+1})\leq G(u^{t+1},u^{t})\leq
G(u^{t},u^{t})\leq F(u^{t})$. Hence, $F$ is nonincreasing under certain update
rules and the convergence of loss (3) is proved.

| Method | ME-LA | PE-BAY | Cov-ca | Cov-de | Movie
---|---|---|---|---|---|---
| G+$D_{I}$+loss 3 | 0.1654 | 0.0793 | 0.0529 | 0.0798 | –
Block-wise | G+$D_{II}$+loss 3 | 0.1708 | 0.0825 | 0.0548 | 0.0879 | –
Missing | G+$D_{I}$+$D_{II}$ | 0.2651 | 0.2082 | 0.1945 | 0.2261 | –
Data | G+$D_{I}$+$D_{II}$+MSE loss | 0.1734 | 0.0818 | 0.0513 | 0.0797 | –
|  Ours(G+$D_{I}$+$D_{II}$+loss 3) | 0.1561 | 0.0743 | 0.0495 | 0.0784 | –
| G+$D_{I}$+loss 3 | 0.1409 | 0.0631 | 0.0346 | 0.0510 | 0.1916
Scattered | G+$D_{II}$+loss 3 | 0.1419 | 0.0693 | 0.0408 | 0.0585 | 0.1921
Missing | G+$D_{I}$+$D_{II}$ | 0.2071 | 0.1602 | 0.1322 | 0.1766 | 0.2810
Data | G+$D_{I}$+$D_{II}$+MSE loss | 0.1398 | 0.0668 | 0.0363 | 0.0564 | 0.1893
|  Ours(G+$D_{I}$+$D_{II}$+loss 3) | 0.1322 | 0.0608 | 0.0321 | 0.0488 | 0.1797
Table 3: RMSE performance in the ablation study. The design of experiments is
the same as 6.1

##  6 Experiments and Evaluation

In this section, we rigorously test the performance of our data imputation
method, BlockEcho, through a series of four carefully designed experiments
using various real-world datasets available in the public domain. The initial
set of experiments compare BlockEcho with state-of-the-art baselines, focusing
on a fixed missing rate as high as 60%. Following this, further experiments
are conducted to assess the stability of these models as the missing rate is
incrementally increased. Additionally, an ablation study is performed to
examine the impact of specific components of the method. Finally, we designed
a downstream prediction task using the imputed data as input, to evaluate how
the performance of the data imputation method influences the end-to-end
performance of a real-world prediction task. The experiments in this project
are based on the following open-source real-world datasets from three distinct
fields:

  1. 1.

_Traffic Datasets (Time-Street matrix)_. We employ two traffic flow datasets,
ME-LA and PE-BAY, which collected in the highway of Los Angeles County and
from California Transportation Agencies Performance Measurement SystemLi et
al. (2017). The traffic data are characterised by strong periodicity such as
peak hours.

  2. 2.

_COVID-19 Dataset (Date-City matrix)_. This government dataset records daily
COVID-19 cases (Cov-ca) and deaths (Cov-de) in major cities around the world
since 2020. Infectious disease datasets are characterized by local continuity
and sudden bursts.

  3. 3.

_Movie Dataset (User-Movie matrix)_. This dataset (Movie) records movie
ratings collected from the MovieLens websiteHarper and Konstan (2015), which
contains a large amount of missing data.

For “scattered missing data,” we generate a random matrix of the same size as
the original data and rank the random numbers within it. We compare the
rankings with $missrate*m*n$ to determine which original elements should be
masked. As for the “block-wise” missing data, we randomly select the top-left
and bottom-right corners of a “block” within the feasible region limited by
the missing rate, and then mask the original elements within that region.

All models were trained on an Nvidia Tesla V100S PCIE GPU and each experiment
is repeated for ten times with different random seeds, and the results are
averaged.

###  6.1 Performance of BlockEcho

In the first series of experiments, we synthetically mask 60% of the data in
each dataset to generate a missing data matrix, and we compare the performance
of BlockEcho with 5 baseline matrix completion models - MissForestStekhoven
and Bühlmann (2011), Matrix FactorizationHastie et al. (2015), GAINYoon et al.
(2018), PC-GAIN Wang et al. (2021) and STGANYuan et al. (2022) — in terms of
imputation accuracy. Data are masked in two ways: block-wise mask and
scattered mask. To measure the imputation accuracy, we use RMSE to calculate
the error between the model-imputed data and the original masked data, which
is defined as
$RMSE=\frac{||(\bar{X}-X)\odot(1-M)||_{F}}{\sum_{i,j}(1-M)_{ij}}$ . In order
to have a relatively uniform measure for different datasets, we normalize the
data before computing RMSE. Table 2 reports the RMSE for BlockEcho and 5 other
imputation models. Results show that BlockEcho performs significantly better
than baselines, especially for block-wise missing data.

###  6.2 Models performance in different missing rates

In the second series of experiments, we take the traffic dataset PE-BAY as a
representative, dynamically adjust the data missing rate from 20% to 80%, and
explore the variation trend of each model performance with the data missing
rate. Figure 3 quantitatively shows this trend: the SOTA machine learning
algorithms represented by GAIN perform better on datasets with a low data
missing rate, but will deteriorate rapidly as the missing rate increases; the
Matrix Factorization method performs more stable but has a lower ceiling for
accuracy. BlockEcho absorbs the advantages of both, and thus gives the best
performance at high missing rates. Although at low missing rates it is not as
accurate as some SOTA models, BlockEcho significantly outperforms them at high
missing rates when the datasets become more incomplete. This superiority is
more obvious in block-wise missing datasets.

###  6.3 Ablation study

In this subsection we design ablation experiments to verify the contribution
of each component of BlockEcho to the results. The table 3 shows that when we
remove any part or replace it with a more conventional alternative, it will
negatively affect the results. It is worth mentioning that when the matrix
factorization loss function 3 is removed and only the GAN framework is used
for matrix completion, the accuracy of the imputed data will be significantly
reduced, which is why most generative models (including ours) use direct loss
functions to guide the convergence of GAN.

###  6.4 Case study: traffic forecasting

In the final series of experiments, we design a downstream prediction task to
illustrate how the performance of the data imputation method affects the end-
to-end performance of a real-world prediction task. We take PE-BAY dataset as
input to forecast the traffic conditions at the next timestamp. For the input
datasets, we take the original data, the imputed data with our model
BlockEcho, and the imputed data from each baseline model for comparison. For
traffic forecasting, we use Random Forests (RF) with the same set of
hyperparameter settings with the same feature engineering for all input
datasets. We use Weighted Mean Absolute Percentage Error (WMAPE) as the error
metric to measure the prediction performance. Table 4 summarizes the
prediction results. Comparing Table 2 and Table 4, we find that models with
higher data imputation accuracy tend to give better results in subsequent
prediction tasks, and BlockEcho, the best-performing data imputation model,
also gives the lowest prediction errors in subsequent traffic forecasting
tasks.

| Block-wise | Scattered
---|---|---
Ori-Data | 0.0378 | 0.0378
MissForest | 0.0696 | 0.0512
Matrix Factorization | 0.0615 | 0.0619
GAIN | 0.0789 | 0.0484
PC-GAIN | 0.0702 | 0.0479
STGAN | 0.0685 | 0.0503
BlockEcho(ours) | 0.0484 | 0.0473
Table 4: WMAPE of prediction task after various imputation model.

##  7 Conclusion and Future Work

In this paper, we excavate and mathematically define the issue of “block-wise”
missing data and innovatively propose the solution, BlockEcho. Experiments on
various data sets, especially with block missing data, show that our method
outperforms other SOTA methods. Our future work will extend to federated
learning where block-wise missing data widely appear.

## References

  * Buuren and Oudshoorn [2000] S. Van Buuren and C. G. M Oudshoorn.  Multivariate imputation by chained equations : Mice v1.0 user’s manual.  2000\.
  * Chang et al. [2023] Zhuoqing Chang, Shubo Liu, Zhaohui Cai, and Guoqing Tu.  Anode-gan: Incomplete time series imputation by augmented neural ode-based generative adversarial networks.  In International Conference on Artificial Neural Networks, pages 16–27. Springer, 2023.
  * Gao et al. [2020] Xin Gao, Fang Deng, and Xianghu Yue.  Data augmentation in fault diagnosis based on the wasserstein generative adversarial network with gradient penalty.  Neurocomputing, 396:487–494, 2020.
  * Gold and Bentler [2000] Michael Steven Gold and Peter M Bentler.  Treatments of missing data: A monte carlo comparison of rbhdi, iterative stochastic regression imputation, and expectation-maximization.  Structural Equation Modeling, 7(3):319–355, 2000.
  * Goodfellow et al. [2020] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial networks.  Communications of the ACM, 63(11):139–144, 2020.
  * Harper and Konstan [2015] F Maxwell Harper and Joseph A Konstan.  The movielens datasets: History and context.  Acm transactions on interactive intelligent systems (tiis), 5(4):1–19, 2015.
  * Hastie et al. [2015] Trevor Hastie, Rahul Mazumder, Jason D Lee, and Reza Zadeh.  Matrix completion and low-rank svd via fast alternating least squares.  The Journal of Machine Learning Research, 16(1):3367–3402, 2015.
  * Heitjan and Basu [1996] Daniel F Heitjan and Srabashi Basu.  Distinguishing “missing at random” and “missing completely at random”.  The American Statistician, 50(3):207–213, 1996.
  * Kong et al. [2023] Xiangjie Kong, Wenfeng Zhou, Guojiang Shen, Wenyi Zhang, Nali Liu, and Yao Yang.  Dynamic graph convolutional recurrent imputation network for spatiotemporal traffic missing data.  Knowledge-Based Systems, 261:110188, 2023.
  * Lee and Seung [2000] Daniel Lee and H Sebastian Seung.  Algorithms for non-negative matrix factorization.  Advances in neural information processing systems, 13, 2000.
  * Lee et al. [2019] Dongwook Lee, Junyoung Kim, Won-Jin Moon, and Jong Chul Ye.  Collagan: Collaborative gan for missing image data imputation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2487–2496, 2019.
  * Li et al. [2017] Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu.  Diffusion convolutional recurrent neural network: Data-driven traffic forecasting.  arXiv preprint arXiv:1707.01926, 2017.
  * Li et al. [2019] Steven Cheng-Xian Li, Bo Jiang, and Benjamin Marlin.  Misgan: Learning from incomplete data with generative adversarial networks.  arXiv preprint arXiv:1902.09599, 2019.
  * Mescheder et al. [2017] Lars Mescheder, Sebastian Nowozin, and Andreas Geiger.  The numerics of gans.  Advances in neural information processing systems, 30, 2017.
  * Miyato et al. [2018] Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida.  Spectral normalization for generative adversarial networks.  arXiv preprint arXiv:1802.05957, 2018.
  * Neves et al. [2021] Diogo Telmo Neves, Marcel Ganesh Naik, and Alberto Proença.  Sgain, wsgain-cp and wsgain-gp: Novel gan methods for missing data imputation.  In International Conference on Computational Science, pages 98–113. Springer, 2021.
  * Santos et al. [2019] Miriam Seoane Santos, Ricardo Cardoso Pereira, Adriana Fonseca Costa, Jastin Pompeu Soares, João Santos, and Pedro Henriques Abreu.  Generating synthetic missing data: A review by missing mechanism.  IEEE Access, 7:11651–11667, 2019.
  * Scheffer [2002] Judi Scheffer.  Dealing with missing data.  2002\.
  * Stekhoven and Bühlmann [2011] D. J. Stekhoven and P Bühlmann.  Missforest - nonparametric missing value imputation for mixed-type data.  2011\.
  * Wang et al. [2021] Yufeng Wang, Dan Li, Xiang Li, and Min Yang.  Pc-gain: Pseudo-label conditional generative adversarial imputation networks for incomplete data.  Neural Networks, 141:395–403, 2021.
  * Williams [2015] R Williams.  Missing data part 1: Overview, traditional methods. university of notre dame, 2015.
  * Yoon and Sull [2020] Seongwook Yoon and Sanghoon Sull.  Gamin: Generative adversarial multiple imputation network for highly missing data.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 8456–8464, 2020.
  * Yoon et al. [2018] Jinsung Yoon, James Jordon, and Mihaela Schaar.  Gain: Missing data imputation using generative adversarial nets.  In International conference on machine learning, pages 5689–5698. PMLR, 2018.
  * Yu et al. [2016] Hsiang-Fu Yu, Nikhil Rao, and Inderjit S Dhillon.  Temporal regularized matrix factorization for high-dimensional time series prediction.  Advances in neural information processing systems, 29, 2016.
  * Yuan et al. [2022] Ye Yuan, Yong Zhang, Boyue Wang, Yuan Peng, Yongli Hu, and Baocai Yin.  Stgan: Spatio-temporal generative adversarial network for traffic data imputation.  IEEE Transactions on Big Data, 2022.
  * Zhang et al. [2021] Weibin Zhang, Pulin Zhang, Yinghao Yu, Xiying Li, Salvatore Antonio Biancardo, and Junyi Zhang.  Missing data repairs for traffic flow with self-attention generative adversarial imputation net.  IEEE Transactions on Intelligent Transportation Systems, 23(7):7919–7930, 2021.


# Heuristic Algorithm for Generalized Function Matching

Radu Stefan Mincu  Department of Computer Science, University of Bucharest,
Bucharest, Romania

###### Abstract

The problem of generalized function matching can be defined as follows: given
a pattern $p=p_{1}\cdots p_{m}$ and a text $t=t_{1}\cdots t_{n}$, find a
mapping $f:\Sigma_{p}\rightarrow\Sigma_{t}^{*}$ and all text locations $i$
such that $f(p_{1})f(p_{2})\cdots f(p_{m})=t_{i}\cdots t_{j}$, a substring of
$t$.

By modifying the restrictions of the matching function $f$, one can obtain
different matching problems, many of which have important applications. When
$f:\Sigma_{p}\rightarrow\Sigma_{t}$ we are faced with problems found in the
well-established field of combinatorial pattern matching. If the single
character constraint is lifted and $f:\Sigma_{p}\rightarrow\Sigma_{t}^{*}$, we
obtain generalized function matching as introduced by Amir and Nor (JDA 2007).
If we further constrain $f$ to be injective, then we arrive at generalized
parametrized matching as defined by Clifford et al. (SPIRE 2009).

There are a number of important applications for pattern matching in
computational biology, text editors and data compression, to name a few.
Therefore, many efficient algorithms have been developed for a wide variety of
specific problems including finding tandem repeats in DNA sequences,
optimizing embedded systems by reusing code etc.

In this work we present a heuristic algorithm illustrating a practical
approach to tackling a variant of generalized function matching where
$f:\Sigma_{p}\rightarrow\Sigma_{t}^{+}$ and demonstrate its performance on
human-produced text as well as random strings.

###### keywords:

string algorithms , pattern matching , heuristics

<EMAIL_ADDRESS>
##  1 Introduction

Given a pattern $p=p_{1}p_{2}\cdots p_{m}$ over an alphabet $\Sigma_{p}$ and a
text $t=t_{1}t_{2}\cdots t_{n}$ over an alphabet $\Sigma_{t}$, the focus of
generalized function matching (GFM) is to find all the occurrences within the
text of the image of the pattern through a function $f$. In other words, the
goal of GFM is to find a function $f:\Sigma_{p}\rightarrow\Sigma_{t}^{*}$ and
all text locations $i$ such that $f(p_{1})f(p_{2})\cdots f(p_{m})=t_{i}\cdots
t_{j}$ i.e. $f(p)$ is a substring of $t$.

### Previous and Related Work

In parametrized function matching introduced by Baker [7], the task is to find
a bijection $f:\Sigma_{p}\rightarrow\Sigma_{t}$. The motivation in Baker’s
work was related to software maintenance: finding pieces of code that are
identical with the exception of variable renaming, with applications in
optimizing embedded systems.

Later on, Amir et al. [1] introduce function matching, a similar problem to
parametrized function matching, where $f:\Sigma_{p}\rightarrow\Sigma_{t}$ can
be any function. The applications include optimizing embedded systems,
techniques for computational biology and image processing.

Following that, Amir and Nor [2] propose generalized function matching (GFM)
where the pattern is used to describe more complicated substructures inside a
text, rather than simply 1:1 substring matching. It is here for the first time
that the image of the matching function is a substring of the text, rather
than a single character i.e. $f:\Sigma_{p}\rightarrow\Sigma_{t}^{*}$.

An $\mathcal{NP}$-completeness result is presented for an altered version of
the problem, the so-called DGFM or generalized function matching with don’t
cares, where wildcard symbols may appear in the text or pattern. Text
wildcards are treated as being equal to all characters of $\Sigma_{t}$, while
pattern wildcards may match any text substring (i.e. the wildcard symbols
found in different pattern locations are considered new, distinct symbols of
$\Sigma_{p}$).

Amir and Nor [2] also describe a simple greedy algorithm to solve GFM that is
exponential in the size of the pattern alphabet $|\Sigma_{p}|$. The algorithm
is introduced out of theoretical interest as it shows the problem is solvable
in polynomial time for bounded alphabets.

In [4] Clifford et al. restrict $f$ to be an injective function and introduce
generalized parametrized matching. Additionally, Clifford et al. produce an
$\mathcal{NP}$-completeness proof for GFM when
$f:\Sigma_{p}\rightarrow\Sigma_{t}^{+}$.

Results regarding the approximability of the problems are given by Clifford
and Popa in [3]. For a detailed survey regarding generalized
function/parametrized matching we refer the reader to Ordyniak and Popa [5,
6].

######  Lemma 1.

Whole text matching is equivalent to substring matching. Finding
$f:\Sigma_{p}\rightarrow\Sigma_{t}^{+}$ such that $f(p)=t_{i}\cdots t_{j}$ is
equivalent to finding
$f:\Sigma_{p}\cup\\{\alpha,\beta\\}\rightarrow(\Sigma_{t}\cup\\{\\#,\$\\})^{+}$
such that $f(p^{\prime})=f(\alpha p\beta)=\\#t\$$, where
$\\{\alpha,\beta\\}\cap\Sigma_{p}=\emptyset$ and
$\\{\\#,\$\\}\cap\Sigma_{t}=\emptyset$.

###### Proof.

The newly introduced pattern symbol $\alpha$ can only be matched to the newly
introduced text symbol $\\#$ or to a non-empty prefix of $t$. Similarly,
$\beta$ can only be matched with $\$$ or a non-empty suffix of $t$. The
remainder of $t$ (which is a substring) has to match the pattern $p$ in the
case of a valid match for $p^{\prime}$. ∎

Using the previous result, we may focus on matching the pattern with the
entire text instead of every text substring.

We make a note here that for the algorithm to be of practical use, we restrict
$f:\Sigma_{p}\rightarrow\Sigma_{t}^{*}$ to
$f:\Sigma_{p}\rightarrow\Sigma_{t}^{+}$. Should we allow empty substrings in
the image of $f$, then we can choose to map variations of $p$ from which we
have removed some symbol(s) from its alphabet (effectively mapping them to the
empty text substring). The practical relevance of this behavior of GFM is
arguable, but we choose to disallow empty substrings to be matched hereafter.

The paper is structured as follows. In Section 2 we give an informal
description for the proposed algorithm and afterwards provide a general
pseudocode implementation and discuss its running time. Section 3 introduces a
heuristic for pattern decomposition that helps reduce the running time when
the pattern has a certain structure. In Section 4 we present two experiments,
one focusing on English sonnets’ structure and the other describing algorithm
running time on random data.

##  2 The Proposed Algorithm for GFM

In this section we detail the inner workings of the heuristic algorithm used
to tackle GFM.

###  2.1 Description of the Algorithm

The input to the algorithm is a finite pattern $p=p_{1}p_{2}\cdots p_{m}$ and
a finite text over some alphabet $t=t_{1}t_{2}\cdots t_{n}$. The output of the
algorithm is a matching function $f:\Sigma_{p}\rightarrow\Sigma_{t}^{+}$ such
that $f(p_{1})f(p_{2})\cdots f(p_{m})=t$.

The key idea of the algorithm is to focus on the repetitions in the pattern
and on the repetitive structure of the input text. The repetitive structure of
the input text can be efficiently discerned by using a suffix tree and this is
the advantage we choose to exploit when designing the present algorithm. What
remains is to understand how we may match the repetitive structure of the
pattern with the repetitive structure of the text.

First, assume that the given pattern has no repetitions. If this is the case
then the problem becomes an exercise in listing all possible factorizations
into $|p|$ non-empty pieces of all substrings of the text. While not entirely
trivial, this is of little practical interest.

Otherwise, there is at least one symbol in the pattern occurring more than
once. In this case, we separate the repeating symbols from the non-repeating
ones. We will refer to the resulting pair of subsequences as the repetitive
subsequence, denoted by $p^{rep}$, and the non-repetitive subsequence, denoted
by $p^{nrep}$, respectively (e.g. if $p=$fabbcdaae, $p^{rep}=$abbaa,
$p^{nrep}=$fcde).

By using the aforementioned repetitive subsequence $p^{rep}$ we can divide the
problem into firstly matching $p^{rep}$ with the text and secondly ensuring
that there exists a valid mapping for the rest of the symbols in the non-
repetitive subsequence $p^{nrep}$. This will include checking whether the
symbols in $p^{nrep}$ are matched to non-empty substrings of the text.

The overall strategy to tackle GFM can now be formulated as follows:

  1. 1.

Discern the repetitive structure of the input text (here we use a suffix
tree).

  2. 2.

Process the pattern, separating a repetitive subsequence and a non-repetitive
subsequence.

  3. 3.

Find the matches of the repetitive subsequence with the text.

  4. 4.

Ensure that the symbols in the non-repetitive subsequence can be matched with
non-empty strings.

###  2.2 The Subsequence Matching in Detail

#### Text Processing

First of all, it is necessary to build a suffix tree for the input text and
store the repetition information into a suitable structure. For the suffix
tree construction step we use an implementation of Ukkonen’s algorithm [8]. To
record a repetition of a substring of the text we have considered the
following structure $(length,occ,[i_{1},\dots,i_{occ}])$, where $length$ is
the number of characters in the repeating substring, followed by a list of
positions where the string appears in the text ($i_{t}$, $2\leq t\leq occ$),
$occ$ being the number of occurrences. In the rest of the paper, the term
repetition will be interpreted as this previously defined structure.

A list of repetition tuples as described above may be obtained by traversing
the constructed suffix tree. Upon traversal, we may choose to discard some
redundant information: the suffixes of a repeating substring that have the
same number of occurrences as the substring. The benefits of this choice
include reducing the search space and memory required for the algorithm.
However, this leads to missing some solutions. To ensure algorithm
correctness, the entire list of repetition tuples is required. In our
heuristic approach, we may choose to search a reduced repetition list for the
first iterations of the algorithm (leading to a much faster search) and then
gradually recover the lost information in subsequent iterations, should no
match be found.

After the repetition list is sorted by length we merge all of the position
lists of the repetitions (and keep a record of the index of the corresponding
repetition and the index of the position inside the position list). In case of
same starting positions, the position of the shorter repeated string appears
first for reasons soon to be discussed in the matching algorithm. This merged
list will be used to match the repetition pattern to the text.

#### Pattern Processing

For the pattern processing step, we first count the occurrences of each
pattern symbol. Symbols that appear multiple times will form the repetitive
subsequence. Additionally, we record the length of each substring made of non-
repeating symbols appearing between each repeating symbol. For an example
consider $p=$fabbcdaae: the repetitive subsequence is $p^{rep}=$abbaa and the
non-repeating string length information can be represented as $count^{nrep}=$
1:0:0:2:0:1.

#### Subsequence Matching

To find a match between the repetitive subsequence and the text, we will map
(substitute) each symbol $\sigma$ in the repetitive subsequence of the pattern
with the index of a repetition having at least as many occurrences as the
occurrences of the symbol in the pattern ($occ\geq|p|_{\sigma}$).

The key observation is that if the repetitive pattern matches the text, then,
for some repetition index substitution, the substituted pattern will be a
subsequence inside the merged repetition list. While this is bound to happen,
it may still be difficult to obtain a valid matching $f$.

#### String Trimming

Once an occurrence of the repetitive subsequence has been detected in the
merged repetition list, it is necessary to “trim” any overlapping substrings.
In short, the starting point of any substring will cut the length of the
previous substring. This change will affect all substrings that belong to a
repetition (they must be identical). If the trimming operation ends up
producing empty strings, then the match is rejected at this point.

Another detail of the trimming process is that it is necessary to sometimes
shorten a repeating string (cut away a suffix) to make room for the non-
repetitive pattern substring in-between. In other words, we must allow the
length of the text substring between repeating pattern symbol matches to be at
least the length of the non-repetitive pattern substring in-between
(effectively ensuring that non-repeating pattern symbols can each be matched
to a non-empty string).

After trimming, we use the information regarding the non-repetitive string
lengths to enforce adjacency between repeating strings should the length be 0.
If this is not possible, reject at this step.

#### Repetition Splitting

If the algorithm reports no match (or insufficient matches) after the previous
steps then it is necessary to reintroduce some of the redundant information
that was discarded during the repetition list construction. To do this, we
simply split each repetition into “halves” by the following method:
$\\{$abcd$\\}$ $\rightarrow$ $\\{$abcd,cd$\\}$ $\rightarrow$
$\\{$abcd,bcd,cd,d$\\}$. Correctness is only ensured by employing the
previously described trimming procedure.

By performing this operation multiple times on the repetition list we can
recover all of the discarded information. The reason for recovering this
information gradually with this heuristic is that the size of the repetition
list greatly impacts the expected running time.

###  2.3 Algorithm pseudo-code

For an insight into the intricacy of the algorithm described in this section,
we refer the reader to the algorithm blocks that contain pseudocode for our
proposed algorithm (Proposed Algorithm block) and for the previous algorithm
by Amir and Nor [2], augmented with some optimizations (Previous Algorithm
block).

Proposed Algorithm: whole-text matches input pattern $p$ with input text $t$
and outputs the respective partitioning of $t$. The notation $|w|_{a}$ denotes
a procedure that returns the number of occurrences of symbol $a$ in string
$w$. The operator $\cdot$ denotes string concatenation. $L^{X}$ is the set of
functions $map:X\rightarrow L$. To keep a high-level view, the complete
trimming procedure is not included (it is lengthy, contains subcases for the
beginning and ending of the text etc.). The variable $subseq$ denotes a list
of positions in the text where the respective subsequence occurs.

input: string $p$,$t$; output: integer $u[1..|p|$]; (such that $t_{u[i]}\cdots
t_{u[i+1]}=f(p_{i})$)

1:$t\leftarrow\\#t\$$, $p^{rep}\leftarrow\emptyset$,
$count^{nrep}_{1}\leftarrow 0$, $j\leftarrow 1$

2:build a suffix tree $ST(t)$

3:traverse $ST(t)$ and build a list $L$ of repetition structures
$(length,occ,pos[1..occ])$

4:$repStruct\leftarrow$ merge and sort the $pos$ arrays of all repetitions in
$L$

5:for $i\leftarrow 1,|p|$ do

6: if $|p|_{p_{i}}>1$ then

7: $p^{rep}\leftarrow p^{rep}\cdot p_{i}$

8: $j\leftarrow j+1$, $count^{nrep}_{j}\leftarrow 0$

9: else

10: $count^{nrep}_{j}\leftarrow count^{nrep}_{j}+1$

11:$X=\Sigma_{p}^{rep}\leftarrow\Sigma(p^{rep})$,

12:for all mapping functions $map\in L^{X}$ do

13: for all occurrences $subseq$ of $map(p^{rep})$ as a subsequence in
$repStruct$ do

14: $L^{copy}\leftarrow L$

15: for $j\leftarrow 1,|subseq|$ do

16: trim the $length$ of $L^{copy}_{i}$ by $subseq_{j+1}-count^{nrep}_{j}$
where $map(p^{rep}_{j})=L^{copy}_{i}$

17: insert the starting positions of the non-repetitive text substrings into
$subseq$

18: if trimming produces no empty substring then

19: output $subseq$

Previous Algorithm: matches input pattern $p$ with input text $t$ and outputs
the respective partitioning of $t$.

input: string $p$,$t$; output: integer $u[1..|p|$]; (such that $t_{u[i]}\cdots
t_{u[i+1]}=f(p_{i})$)

1:$i\leftarrow 1$

2:for $c\in\Sigma_{p}$ do

3: $index_{c}\leftarrow i$, $lengths_{i}\leftarrow 1$, $i\leftarrow i+1$

4:for $i\leftarrow 1,|t|$ do

5: $availableLen\leftarrow|t|-i$

6: for every combination of $lengths$ s.t. $sum(lengths)\leq availableLen$ do

7: for $j\leftarrow 1,|\Sigma_{p}|$ do

8: $starts_{j}\leftarrow=-1$

9: $OK\leftarrow true$, $last\leftarrow i$

10: for $j\leftarrow 1,|p|$ do

11: $c\leftarrow p_{j}$, $k\leftarrow index_{c}$

12: if $starts_{k}=-1$ then

13: $starts_{k}\leftarrow last$

14: else

15: for $\ell\leftarrow last,last+lengths_{k}$ do

16: if $t_{\ell}\neq t_{\ell+start_{k}-last}$ then

17: $OK\leftarrow$ false, break

18: if $OK=$ false then

19: break

20: $last\leftarrow last+lengths_{k}$

21: $last=i$, output $last$

22: for $j\leftarrow 1,|p|$ do

23: $c\leftarrow p_{j}$, $k\leftarrow index_{c}$

24: $last\leftarrow last+lengths_{k}$, output $last$

###  2.4 Time complexity analysis

In the upcoming analysis we consider $|t|=n$ and $|p|=m$. Additionally we use
the denomination uniformly random string or uniform string when we refer to a
string that is built by concatenating symbols in a finite, given alphabet and
each symbol has equal probability to be chosen for each position in the
string.

From the proposed algorithm’s pseudocode we observe that the computation time
is dominated by the two loops in lines 11-12. Before that, a suffix tree is
built in linear time $O(n)$ using Ukkonen’s algorithm, then the list of
repetitions $L$ is recovered by traversing the tree in $O(n)$. Following, the
repetition list has $O(n)$ structures each containing $O(n)$ integers in their
list of occurrences (in the case of uniformly random strings the expected
length of these lists is $O(n|\Sigma_{t}|^{-1})$). As a result, merging and
sorting the lists of occurrences takes up to $O(n^{2}\log n)$ time (and
$O(n^{2}|\Sigma_{t}|^{-1}\log n^{2}|\Sigma_{t}|^{-1})$ expected time for
uniform strings) using comparison based sorting. The expected size of the
sorted occurrence list is therefore between $O(n^{2}|\Sigma_{t}|^{-1})$ in the
uniform string case and $O(n^{2})$ in general . The pattern processing step is
completed in $O(m)$.

The dominating term of the time complexity of the proposed algorithm is
depending on the size of the repetition list. The number of
$map:\Sigma^{rep}_{p}\rightarrow L$ functions is $O(n^{|\Sigma_{p}|}$). For
each such function, we search for a subsequence of length $|p^{rep}|$ inside
the merged repetition occurrence list.

The search for the first occurrence of a subsequence in the lexicographic
order can be done in linear time with regard to the size of the occurrence
list (uniform string case $O(n^{2}|\Sigma_{t}|^{-1})$, in general $O(n^{2})$).
The next subsequence ocurrence in the lexicographic order of positions can be
found in time $O(mn^{2}|\Sigma_{t}|^{-1})$ for the uniform string case,
$O(mn^{2})$ in general). In the worst case, there are generally ${n^{2}\choose
m}=O(n^{2m})$ subsequence occurrences and $O(n^{2m}|\Sigma_{t}|^{-m})$ for
uniform strings. For each occurrence a copy of the repetition list is made.
Actually, only the repetitions that appear in the image of the mapping
function need to be duplicated, so this step takes only $O(|\Sigma_{p}|)$.

The trimming of the length of each repetition occurrence against the starting
position of the next occurrence will take constant time, so the entire
trimming procedure takes time $O(m)$, including the handling of the non-
repetitive subsequence.

With the previous explanation, the total time complexity of the algorithm is
generally $O(mn^{2m+|\Sigma_{p}|})$ and for uniform strings
$O(mn^{2m+|\Sigma_{p}|}\cdot|\Sigma_{t}|^{-m})$.

The reference implementation of the algorithm by Amir and Nor has been
optimized to $O(n^{|\Sigma_{p}|+1})$.

##  3 Pattern decomposition heuristics

Our search for matches is based around finding a subsequence within the
repetition structure obtained from all possible substitutions of the
repetitive pattern subsequence with a valid assignment. As such, even for
relatively small tests, if the text has a highly repetitive structure, then
the search is prohibitively slow.

Our proposed approach to mitigate this drawback is a divide-and-conquer
heuristic. The idea is that if a pattern can be split into pieces that have
disjoint alphabets, then we can try to match once of the pieces to obtain a
candidate solution fragment to the whole pattern match problem. Subsequently,
for each candidate fragment identified, we solve the problem of matching the
remaining pieces of the pattern with the remaining pieces of the text.

To illustrate our approach consider the following example. Let
$p=ABABCDCDEFEFGG$ the rhyme structure of Shakespearean sonnets (see
Subsection 4.1 for an experiment regarding this pattern). We can split $p$ in
the following format to facilitate our search: $p=ABAB\cdot CDCD\cdot
EFEF\cdot GG$. In other words, first we locate an occurrence of $ABAB$, then
we find a match of $CDCD$ that begins on the first position of the remaining
text, and so on.

This pattern splitting method is better formalized by using whole-text
matching. If we try to match $p=*ABABCDCDEFEFGG*$ to the entire text, then we
can easily represent the split pattern as $(*ABAB*,CDCD*,EFEF*,GG*)$. The
previously described technique only works left-to-right: first $*ABAB*$ is
matched to the text, then $CDCD*$ is matched to the text substring that
matches the rightmost $*$ symbol from $*ABAB*$, then $EFEF*$ is matched to the
text substring matching $*$ from $CDCD*$, and so on.

So far we have described a type of pattern splitting of $p$ into as many
pieces $q_{i}$ as possible such that $p=q_{1}\cdots q_{k}$ and all $q_{i}$
have disjoint alphabets i.e. $\Sigma(q_{i})\cap\Sigma(q_{j})=\emptyset$,
$\forall i\neq j$. This method is based on splitting the pattern in the manner
$p=qr$ with $\Sigma(q)\cap\Sigma(r)=\emptyset$. A more general version of this
pattern splitting method is to reduce the pattern “inwards” as inspired by
palindrome type patterns. This method will separate $p=qp^{\prime}r$ with
$p^{\prime}\cap(q\cup r)=\emptyset$. In this case, a whole-text match of
$p=*qp^{\prime}r*$ will be separated into matching $*q\\!*\\!r*$ first and
subsequently matching $p^{\prime}$ with the text substring matching the middle
wildcard of $*q\\!*\\!r*$.

##  4 Experimental results

###  4.1 A toy experiment

To showcase the proposed algorithm, we prepare a toy example inspired from
[2], where the authors give music and poetry analysis as a motivation for GFM
and give the example of the Shakespearean sonnet: a poem with a very rigid
structure that respects the rhyme template $ABABCDCDEFEFGG$.

In this first experiment we use as a pattern $p=ABABCDCDEFEFGG$ and as a text
we randomly select 6 sonnets from the work of Shakespeare, concatenate them,
transcribe them to IPA (International Phonetic Alphabet) and keep only the
last character of each verse.

In Figure 1 we use a graphical representation of the text by using a technique
similar to the one proposed by Wattenberg [9] namely the arc diagram, where
the ends of each arc indicate exact occurrences of a substring and the arc
width represents the length of the repeated substring. Below the gray
repetition structure we also display the matches found with our algorithm
using red hues. In blue we show a correct match that is not a sonnet (but has
a structure that matches the pattern).

It is a noteworthy and surprising find that the third sonnet is not a match
for the given pattern. Upon closer inspection the culprit is Shakespeare’s
Sonnet XLIX, which uses an imperfect rhyme in verses 5-8. Thus, the mismatch
is not through any fault of the algorithm (it is working correctly for proper
matches).

Fig. 1: Graphical representation of the repetitive structure (in grey) of a
string obtained by concatenating the last character in the International
Phonetic Alphabet transcription of each verse of 6 concatenated Shakespearean
sonnets (delimited by the red markers). Mirrored and in red hues are the
matches of the Shakespearean sonnet rhyme pattern $ABABCDCDEFEFGG$ as produced
by our proposed algorithm. In blue is an example of a correct pattern match
that is not a sonnet itself, but is an effect of concatenation.

###  4.2 Running time analysis

For the purpose of determining running time, our proposed algorithm against is
tested against our own optimized implementation of the previous Greedy
algorithm proposed by Amir and Nor [2].

The dataset used in the experiments contains uniformly random strings (meaning
that each position in an uniform string can be any symbol in a fixed alphabet
$\Sigma$ with equal probability $\frac{1}{|\Sigma|}$). The test patterns have
the additional restriction that the first and last symbol must repeat in the
pattern (at least once) and there can be no two consecutive non-repeating
symbols. If this were not the case then the number of valid solutions can
needlessly multiply by $O(|t|)$ e.g. when $p=ab$ is matched to any text $t$
then $a$ is matched to $|t|-1$ prefixes of the text.

We summarize the results of the experiments using Figures 2 and 3: on the
$x$-axis are separate test instances, while on the $y$-axis we plot the
recorded wall clock time. Both algorithms are run on every test instance and
we do not employ pattern decomposition heuristics (from Section 3) for a
proper comparison. The results in Figures 2 and 3 are for 500 instances each.

In Figure 2 we observe the exponential growth in time required for the
proposed algorithm with regard to the amount of distinct repeating substrings
in the text (which grows inversely proportional to the size of the text
alphabet). We also showcase the differences induced by different patterns,
which end up being quite significant. Here we analyze the algorithms on all
patterns of length 5 and alphabet size 3 that respect the previously mentioned
restriction (amounting to 10, disregarding isomorphism). The previous
algorithm by Amir and Nor seems mostly unaffected by text repetition, showing
only a weak linear growth with regard to this parameter (which is in the range
60-130).

In Figure 3 we show the exponential growth in the time of the algorithm by
Amir and Nor with regard to pattern alphabet size. All of the test instances
use the same text and the pattern alphabet varies. We sort the results by
pattern alphabet size and execution time of the previous algorithm. While the
greedy algorithm is quite predictable and is highly dependent on pattern
alphabet, there is much greater variance in the time required for the proposed
algorithm, which is highly dependent on the pattern structure as well
(adjacency of repeating pattern symbols etc.). For all of the instances, we
stop the algorithms when the time exceeds 1 second. It is noteworthy that the
expected running time for the algorithm of Amir and Nor for pattern alphabet
size 4 is around 21-25 seconds.

It is interesting that if we keep the best time of the two algorithms, the
test instances are all solved under 1 second.

Fig. 2: Algorithm running time for uniformly random pattern strings and
uniformly random text strings with the following parameters: fixed pattern
length $|p|=5$, fixed pattern alphabet size $|\Sigma_{p}|=3$, fixed text
length $|t|=200$, variable text alphabet size
$|\Sigma_{t}|\in\\{4,\dots,16\\}$. Fig. 3: Algorithm running time for
uniformly random pattern strings and only one uniformly random text string
with the following parameters: fixed text (fixed text length $|t|=200$, fixed
text alphabet size $|\Sigma_{t}|=16$, distinct repeating substring count is
70), fixed pattern length $|p|=6$, variable pattern alphabet size
$|\Sigma_{p}|\in\\{2,3,4\\}$.

##  5 Conclusion and open problems

We show that while the algorithm by Amir and Nor works well for patterns with
a small alphabet, there are cases where our algorithm is better, notably on
text that has a low number of repeating substrings (which is commonly the case
for English text, since the alphabet is not small enough to prohibit our
proposed algorithm). We also notice that there are types of patterns for which
our algorithm outperforms the previously known algorithm.

For future work it is of interest to design more algorithms for generalized
function matching and develop a specialized toolkit that takes into account
the pattern structure and the text repetitiveness and automates the choice of
the appropriate algorithm.

## References

  * Amir et al. [2003] Amir, A., Aumann, Y., Cole, R., Lewenstein, M., Porat, E., 2003.  Function matching: Algorithms, applications, and a lower bound, pp. 929–942.  doi:10.1007/3-540-45061-0\\_72.
  * Amir and Nor [2004] Amir, A., Nor, I., 2004.  Generalized function matching, pp. 361–378.  doi:10.1007/978-3-540-30551-4\\_6.
  * Clifford and Popa [2010] Clifford, R., Popa, A., 2010\.  (in)approximability results for pattern matching problems, pp. 52–62.
  * Clifford et al. [2009] Clifford, R., Wettroth Harrow, A., Popa, A., Sach, B., 2009.  Generalised matching, pp. 295–301.  doi:10.1007/978-3-642-03784-9\\_29.
  * Ordyniak and Popa [2014] Ordyniak, S., Popa, A., 2014\.  A parameterized study of maximum generalized pattern matching problems, in: Parameterized and Exact Computation, Springer International Publishing, Cham. pp. 270–281.  doi:10.1007/978-3-319-13524-3\\_23.
  * Ordyniak and Popa [2015] Ordyniak, S., Popa, A., 2015\.  A parameterized study of maximum generalized pattern matching problems.  Algorithmica 75, 1–26.  doi:10.1007/s00453-015-0008-8.
  * S. Baker [1993] S. Baker, B., 1993.  Theory of parameterized pattern matching: Algorithms and applications, pp. 71–80.  doi:10.1145/167088.167115.
  * Ukkonen [1995] Ukkonen, E., 1995.  On-line construction of suffix trees.  Algorithmica 14, 11\.  doi:10.1007/BF01206331.
  * Wattenberg [2002] Wattenberg, M., 2002.  Arc diagrams: Visualizing structure in strings, pp. 110 – 116.  doi:10.1109/INFVIS.2002.1173155.