# Spatially and Robustly Hybrid Mixture Regression Model for Inference of
Spatial Dependence

Wennan Chang1, Pengdao Dang1, Changlin Wan1, Xiaoyu Lu2, Yue Fang3,
Tong Zhao4, Yong Zang3, Bo Li5, Chi Zhang1 and Sha Cao3 1Electrical and
Computer Engineering, Purdue University 2School of Informatics and Computing,
Indiana University
3Department of Biostatistics, Indiana University, Indianapolis, USA, 46202
4Amazon, Seattle, USA 5School of Economics, Peking University, Beijing, China
Corresponding Author Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a Spatial Robust Mixture Regression model to
investigate the relationship between a response variable and a set of
explanatory variables over the spatial domain, assuming that the relationships
may exhibit complex spatially dynamic patterns that cannot be captured by
constant regression coefficients. Our method integrates the robust finite
mixture Gaussian regression model with spatial constraints, to simultaneously
handle the spatial nonstationarity, local homogeneity, and outlier
contaminations. Compared with existing spatial regression models, our proposed
model assumes the existence a few distinct regression models that are
estimated based on observations that exhibit similar response–predictor
relationships. As such, the proposed model not only accounts for
nonstationarity in the spatial trend, but also clusters observations into a
few distinct and homogenous groups. This provides an advantage on
interpretation with a few stationary sub-processes identified that capture the
predominant relationships between response and predictor variables. Moreover,
the proposed method incorporates robust procedures to handle contaminations
from both regression outliers and spatial outliers. By doing so, we robustly
segment the spatial domain into distinct local regions with similar regression
coefficients, and sporadic locations that are purely outliers. Rigorous
statistical hypothesis testing procedure has been designed to test the
significance of such segmentation. Experimental results on many synthetic and
real-world datasets demonstrate the robustness, accuracy, and effectiveness of
our proposed method, compared with other robust finite mixture regression,
spatial regression and spatial segmentation methods.

###### Index Terms:

robust mixture regression, spatial information, hybrid, Markov random field
(MRF), spatial inference

##  I Introduction

Many problems in the environmental, economic, and biological sciences involve
spatially collected data, and a main problem of interest is investigation of
the relationship between a response variable and a set of explanatory
variables over the spatial domain using regression modeling. Notably, the
relationships between response variables and covariates may exhibit complex
spatially dynamic patterns that cannot be captured by constant regression
coefficients. Instead, such relationships may abruptly change at a certain
boundary of two neighboring spatial clusters, but stay relatively homogeneous
within clusters. Detecting clusters of observations that display similarity in
both regression relationships and spatial proximity allows straightforward
interpretations of local associations between response variables and
covariates. For example, the residential real estate pricing could be quite
similar in a local community, but drastically differ for two houses across the
street [1]; a major goal of analyzing functional magnetic resonance imaging
(fMRI) data is to detect spatially distributed and functionally linked regions
that continuously share information with each other in reaction to different
stimuli [2]. For all these real-world application settings, the collected data
may often contain outliers, which may severely corrupt the analysis results if
not properly handled. Overall, the spatial nonstationarity, local homogeneity,
and model robustness are three main challenges in spatial regression modeling.

In the nonspatial setting, finite mixture regression models have been used in
many areas as an effective exploratory approach to identify heterogeneity in
response–predictor relationships. For an overview, see [3, 4]. To account for
outliers or heavy-tailed noises, many algorithms have been developed to
estimate the parameters robustly [5]. To seek for robust parameter estimation
in the presence of outliers, methods have been developed that replaced the
least-square criterion in the M-step of the expectation maximization (EM)
algorithm by more robust criterion [6, 7, 8, 9, 10, 11, 12]. To enable
simultaneous model estimation and outlier removal, penalized mean-shift
mixture model [13], and the least trimmed likelihood estimator [14, 15, 16,
17] were proposed. While these methods could robustly capture the
heterogeneous relationship between response and predictor variables, they are
not designed to model the spatial dependency.

In modeling the spatial dependency, conventional nonstationary spatial
regression models such as geographically weighted regression (GWR) [18, 19,
20] and Bayesian spatially varying coefficient (SVC) [21, 22] models fit as
many regression models to the data as there are observations, at the cost of a
large computational burden for large spatial datasets, and sometimes may lead
to overfitting. In addition, interpretation of the GWR and SVC models require
visual inspection of the coefficient maps to pursue local homogeneity, and can
not automatically capture the spatially clustered patterns. In order to
automatically detect spatially homogeneity cluster, a penalized spatial
regression model has been proposed [23], where a fused-lasso [24] type of
penalty has been developed to account for the spatial homogeneity in the
linear regression setting. Nevertheless, the spatial smoothness assumption in
the above spatial regression models could be problematic and violated due to
natural or man-made discontinuities in the spatial domain. In addition, none
of these methods is designed to handle outliers.

Model-based spatial segmentation is another type of methods to deal with
spatial data using spatially constrained Gaussian mixture model [25, 26].
Spatial segmentation incorporates spatial information between neighboring
pixels into the Gaussian mixture model based on Markov random field (MRF),
with a goal to cluster all variables (e.g. pixels in image), where the
distance of two instances is dependent on both their feature expressions and
spatial proximity. This comes at a high computational cost. While robust
spatial segmentation algorithms are available [25, 26], they fail to
intentionally model the linear relationship between the response and
predictors, but instead simply treat the response and predictors as different
features.

In summary, none of the existing methods could robustly model the spatial
clustering patterns of linear dependency between response and predictors, and
we propose a novel Spatial Robust Mixture Regression (SRMR) model that enables
a simultaneous detection of spatial regions in which variables have a strong
linear dependency.

The key contributions of work include: (1) We developed the very first
computational concept of spatially dependent mixture regression analysis. (2)
We provided the SRMR model that efficiently solves the spatially dependent
mixture regression problem, which is also empowered by a statistical inference
approach to assess regression significance. (3) SRMR enables a new type of
spatial segmentation analysis to detect overlapped spatial regions of varied
dependencies among subset of features, which have high contextual
meaningfulness.

##  II Preliminary

###  II-A Notations

We denote scalar value, vector, and matrix as lowercase character $x$, bold
lowercase character x, and uppercase character $X$, respectively. Let
$\\{(\boldsymbol{x}(s_{i}),y(s_{i})),i=1,\ldots,n\\}$ represent a set of
spatial data that is observed at spatial locations
$s_{1},\ldots,s_{n}\in\mathbb{R}^{2}$, where the response variable $y(s_{i})$
is assumed to be spatially correlated, $\boldsymbol{x}(s_{i})=$
$(x_{1}(s_{i}),\ldots,x_{p}(s_{i}))^{\mathrm{T}}$ is the $p$-dimensional
vector of explanatory variables for the observation located at $s_{i}$, and
$s_{i}=(c_{i}^{1},c_{i}^{2})$ is the 2-dimensional coordinate of the $i$th
location. In this study, we only describe and validate the SRMR model on
2-dimensional spatial data. Noted, the approach can be directly applied to
$K$-dimensional ($K>2$) spatial data.

###  II-B Problem statement

To capture the spatially dependent structure for the response variable, we
write a standard generalized linear regression model (GLM) for the $i$-th
spatial location as follows,

$g(E(y=y(s_{i})\mid\boldsymbol{x}=\boldsymbol{x}(s_{i})))=\sum_{j=1}^{p}x_{j}(s_{i})\beta_{ji}+\epsilon_{i}$

where $\beta_{ji},j=1,2,\ldots,p$, are the regression coefficients for the $p$
predictors, and $\epsilon_{i}$ represents random noise with mean 0 and
variance $\sigma_{i}^{2}$, and $g(\cdot)$ is the link function. In this work,
we assume identify link for linear regression. The intercept can be
accommodated by including 1 as an entry of
$\boldsymbol{x}(\boldsymbol{s}_{\boldsymbol{i}})$. Apparently, unless with
sufficient number of repeated measurements for each location, the
$\beta_{ji},\sigma_{i}$ are non-identifiable. In many cases, there is only a
single observation for each spatial location, certain spatial constraints will
be enforced to ensure the identifiablity of the model parameters.

###### Definition 1

Spatially Dependent Mixture Regression. Given a dataset consisting of $n$
observations $\\{(\boldsymbol{x}(s_{i}),y(s_{i})),i=1,\ldots,n\\}$ from
spatial locations $s_{1},\ldots,s_{n}$, the goal of spatially dependent
mixture regression is to identify spatial regions $\Pi_{1},\ldots,\Pi_{K}$ and
the number $K$, s.t.,

$y(s_{i})=\sum_{j=1}^{p}x_{j}(s_{i})\beta_{j}^{k}+\epsilon_{i},\ if\
s_{i}\in\Pi_{k}$

, where $\beta_{j}^{k},j=1,...,p,k=1,...,K$ are regression parameters for the
$p$ predictors in the $k$-th cluster;
$\epsilon_{i}\sim\mathcal{N}(0,\sigma_{k})$, where $\sigma_{k}$ represents the
noise level of cluster $k$.

To account for the presence of outliers, we assume that $\Pi_{K}$ are non-
overlapping subsets of the whole set $\\{1,...,n\\}$, and denote the outlier
set as $\Pi_{0}$, such that
$\Pi_{0}=\\{1,...,n\\}\setminus\bigcup_{k=1}^{K}\Pi_{k}$ Two type of outliers
will be considered here:
Type 1 Outliers:
$y(s_{i})\not=\sum_{j=1}^{p}x_{k}(s_{i})\beta_{j}^{k}+\epsilon^{k},\forall
k=1,...,K$
Type 2 Outliers: $\exists
k,y(s_{i})=\sum_{j=1}^{p}x_{k}(s_{i})\beta_{j}^{k}+\epsilon^{k},s_{i}\not\in\Pi_{k}$
Here the Type 1 Outliers represent the samples do not fit any regression model
while the Type 2 Outliers represent the ones fit a certain model but do not
locate nearby the spatial region.

Noted, pre-assumptions of the spatial regions $\Pi_{1},\ldots,\Pi_{K}$ are
needed to enable a valid solution of the spatially dependent mixture
regression problem. Such assumptions include a connected spatial region, a
compact shape, or high enrichment to a certain region. Noted, as spatially
dependent mixture regression assigns each sample into one spatial region
$\Pi_{k}$, it directly forms a spatial segmentation method.

###  II-C Related works

Mixture regression and robust estimators. Consider an finite mixture Gaussian
regression model parameterized by
$\boldsymbol{\theta}=\\{(\pi_{k},\boldsymbol{\beta}_{k},\sigma_{k}^{2})\\}_{k=1}^{K}$,
the conditional density of $y$ given $\boldsymbol{x}$ is
$f(y|\boldsymbol{x,\theta})=\sum_{k=1}^{K}\pi_{k}\mathcal{N}(y;\boldsymbol{x^{T}\beta_{k}},\sigma_{k}^{2})$,
where $\mathcal{N}(y;\mu,\sigma^{2})$ is the normal density function with mean
$\mu$ and variance $\sigma^{2}$. Many algorithms have been developed to
estimate the parameters robustly [5] by replacing the least-square criterion
in the M-step with more robust criterion in the EM algorithm [8, 10, 11, 12].
To enable simultaneous model estimation and outlier detection, usually a
hyperparameter regarding the proportion of outlying samples needs to be
specified, such as in the penalized mean-shift mixture model [13], and the
least trimmed likelihood estimator [14, 27].

Spatially smooth regression. Conventional nonstationary spatial regression
models such as geographically weighted regression (GWR) [18, 19, 20] and
Bayesian spatially varying coefficient (SVC) [21, 22] models allow regression
coefficients to vary smoothly as a function of the spatial domain. For GWR,
assuming a linear model with $\boldsymbol{y}$ denoting the observed response
vector and $\boldsymbol{X}$ the design matrix, the regression coefficient at
the $i$ th location is estimated from
$\widehat{\boldsymbol{\beta}}_{i}=(\boldsymbol{X}^{\mathrm{T}}\boldsymbol{W}_{i}\boldsymbol{X})^{-1}\boldsymbol{X}^{\mathrm{T}}\boldsymbol{W}_{i}\boldsymbol{y}$,
where $\boldsymbol{W}_{i}$ is a diagonal weight matrix defined by a kernel
function of distance of all other points to point $i$. The challenge with GWR
and SVC models is that they fit as many regression models to the data as there
are observations, at the cost of a large computational burden, possible over-
fitting and interpretation. A penalized spatial regression model has been
developed to automatically detect clusters [23] by incorporating a fused-lasso
penalty constructed based on spatial proximity.

Spatial segmentation. Model-based spatial segmentation aims perform a
segmentation task to all samples (e.g. pixels in image) based on the input
features. Model-based spatial segmentation adopts an energy function $U(\Pi)$
to integrate the spatial information such as neighborhoods with a regular
clustering analysis of the features. Intrinsically, such methods leverage
spatial and data consistency to segment spatial regions, i.e., only
considering the covariance of independent variables, which cannot solve the
spatially dependent regression problem.

##  III Method

To solve the problem of spatially dependent mixture regression, computational
challenges arise from three aspects: (1) the mixture regression model and
spatial consistency do not form one unified likelihood function, which
prohibits a direct solution by using EM algorithm, (2) detection spatial
regions should depend on both goodness of fitting and spatial consistency, and
(3) there is lack of a validate approach to assess the statistical
significance of mixture regression models.

Input: Response vector ${Y}$; independent variables in matrix
$X^{N\times(P+1)}$; the number of mixing component $K$; size of initialization
random sample $n_{0}$; 2-dimentional spatial coordinates $S^{N\times 2}$;
hyperparameter $\lambda$.

Output: Partition $\mathcal{C^{*}}=\bigcup_{k=1}^{K}C_{k}$; Mixture regression
model parameter estimate $\boldsymbol{\theta}^{*}$; spatial centriod parameter
$w^{*}$ ; type 2 outlier set $U^{*}$

Initialization:
$\mathcal{C}=\\{C_{1},...,C_{K}\\},C_{i}\subseteq\\{1,...,N\\}$ based on
coordinate $S$; compute centroid point $w=\\{w_{1},...,w_{K}\\}$ with
$\mathcal{C}$

1

2for  _$m=0,...,L_{0}$ or until convergence_ do

3

4 E-step: Compute for $i=1,...,N$ and $k=1,...,K$, the hybrid posterior
probabilities $p_{ik}^{(m)}$ by

5
$p_{ik}^{(m)}=(1-\lambda)p_{reg}(z_{i}=k|\boldsymbol{x_{i}},y_{i},\boldsymbol{\theta}^{(m)})+\lambda
p_{spa}(z_{i}=k|S_{i,:},w^{(m)})$

6

7 C-step: For $i=1,...,K$, assign
$C_{k}^{(m)}=\\{i|\underset{l\in\\{1,...,K\\}}{\operatorname{argmax}}\
p_{il}^{(m)}=k,i=1,...,N\\}$, and let $n_{k}^{(m)}$ be the size of
$C_{k}^{(m)}$, $U_{k}^{(m)}=\\{i|z_{ir}\neq
z_{is},z_{ir}=\underset{z_{ir}\in\\{1,...,K\\}}{\operatorname{argmax}}p_{reg},z_{is}=\underset{z_{is}\in\\{1,...,K\\}}{\operatorname{argmax}}p_{spa}\\}$

8

9 M-step: For $k=1,...,K$, the parameters are then updated by
$\pi_{k}^{(m+1)}=\frac{n_{k}^{(m)}}{\sum_{l=1}^{K}n_{l}^{(m)}}$,
$(\boldsymbol{\beta}_{k}^{(m+1)},\sigma_{k}^{2(m+1)})=\textbf{OLS}({{Y}_{C_{k}^{(m)}},X_{C_{k}^{(m)},:}})$
, $w_{k}=\frac{1}{n_{k}^{(m)}}\sum_{i\in C_{k}}S_{i,:}$

10 end for

Algorithm 1 Hybrid Mixture Regression (HMR)

###  III-A SRMR algorithm and mathematical considerations

In sight of the challenge, we developed the spatial robust mixture regression
(SRMR) algorithm to conduct simultaneous outlier detection and spatially
dependent mixture regression estimation. The underlying idea is that by
assuming a likelihood function of spatial regions $p_{spa}$ and introducing a
tuning parameter $\lambda\in(0,1)$ to link $p_{spa}$ with the likelihood of
mixture regression $p_{reg}$, a surrogate likelihood function
$(1-\lambda)p_{reg}+\lambda p_{spa}$ is developed to enable a modified EM-
algorithm (Algorithm 1). The inputs of Algorithm 1 include the response and
independent variables, spatial coordinates, and the hyper parameter $\lambda$.
It conducts a simplified spatially dependent mixture regression fitting by
assuming there is only Type 2 outliers, i.e., the sample fit one mixture model
but do not locate in the corresponding spatial region. Hence, Algorithm 1 fits
a conventional mixture regression model and computes the spatial regions that
are top enriched by the samples fit each regression component. In this study,
we assume the spatial likelihood follows
$p_{spa}(z_{i}=k|s_{i},w)\propto||s_{i},w||_{2}$, where $z_{i}$ represents the
class of sample $i$ and $||s_{i},w||_{2}$ represents the Euclidean distance
between the spatial coordinate of the sample $s_{i}$ and the centers of the
spatial regions $w$, i.e., assuming the spatial regions form a compact shape.
Specifically, a voting step (C-step) is introduced in Algorithm 1, which
identifies Type 2 outliers by the ones whose most likely regression component
and spatial region are not consistent. Noted, as all the input samples are
utilized in the estimation of the mixture regression model, Algorithm 1 is
always convergent.

Based on the Algorithm 1, we developed the SRMR framework (Algorithm 2). In
SRMR, we iteratively conduct the Type 2 outlier only spatially dependent
mixture regression by using the Algorithm 1 and identify Type 1 outliers by
running a robust linear regression on all the samples predicted to each
spatial region. The underlying consideration is that only one regression
component is consisted within each identified spatial region, which could be
effectively identified by a conventional robust regression approach (RLM). In
SRMS, we implement the trimmed likelihood estimation based robust mixture
regression. The inputs of SRMR is the same as the input of Algorithm 1 plus
the maximal iteration number $L_{0}$ and a random seed. The outputs of SRMR
include the identified mixture regression models and outliers. The component
of each non-outlier samples can be further assigned by maximal likelihood. In
SRMR, we utilize the same BIC function for conventional robust mixture
regression analysis.

###  III-B Statistical Inference

####  III-B1 Hypothesis testing for spatial regions

We conducted a geometry based approach to estimate the significance to observe
a spatial region of a certain size. Noted, we utilized the compact spatial
shape assumption in SRMR, which could be considered as a round shape. For a
round shape with a diameter of $r$, the number of the shapes needed to cover a
rectangular spatial region can be computed by $0.28m\times n/r^{2}$, which
serves as a weight to correct the p value assessed from each single component
robust regression as detailed following.

####  III-B2 Hypothesis testing for robust linear regression

We discuss hypothesis testing of the significance a robust linear regression
model parameterized by
$\hat{\boldsymbol{\theta}}=\\{\hat{\boldsymbol{\beta}},\hat{\sigma},\hat{\boldsymbol{\eta}}\\}$,
which represents the robust regression coefficients estimator, standard
deviation estimator, and the index of the outlying samples respectively. A
bootstrap procedure is adopted to test the null hypothesis
$\boldsymbol{\theta}=\hat{\boldsymbol{\theta}}$. We perform the following
steps.

Step 1: Calculating the residuals for all observations, including the outlier
samples, under regression parameter $\hat{\boldsymbol{\beta}},\hat{\sigma}$,
denoted as $\boldsymbol{\epsilon}=\\{\epsilon_{1},...,\epsilon_{n}\\}$. Let
$\boldsymbol{\epsilon}_{out}$ be the residuals corresponding to outlying
samples, and $\epsilon_{0}$ be smallest absolute residual in
$\boldsymbol{\epsilon}_{out}$.

Step 2: Generate iid sample $\tilde{\epsilon}_{1},\ldots,\tilde{\epsilon}_{n}$
from the normal distribution $\mathcal{N}(0,\hat{\sigma})$, denoted as
$\tilde{\boldsymbol{\epsilon}}=(\tilde{\epsilon}_{1},\ldots,\tilde{\epsilon}_{n})$.

Step 3: Calculate the percentage of samples in $\tilde{\boldsymbol{\epsilon}}$
whose absolute values are larger than $\epsilon_{0}$, and denote it as
$p_{0}$.

Step 4: Repeat steps 2-3 for $B$ times, and the statistical significance is
evaluated as the average of $p_{0}$ for the $B$ times.

###  III-C Discussion

Several prominent features make our proposed approach attractive. First
instead of using a robust estimation criterion or complex heavy-tailed
distributions to robustify the mixture regression model, our method is built
upon a spatial regression model so as to facilitate computation and model
interpretation. Second we adopt a sparse and scale-dependent mean-shift
parameterization. Each observation is allowed to have potentially different
outlying effects across different regression components, which is very
flexible. Compared to existing spatial regression methods, our approach allows
an efficient solution via the celebrated penalized regression approach, and
different information criteria (such as AIC and BIC) can be used to adaptively
determine the proportion of outliers. In the next section, we utilized
extensive simulations to demonstrate the performance of SRMR and its highly
robustness to both gross outliers and high leverage points.

Input: Response vector ${Y}$; independent variables in matrix
$X_{N\times(P+1)}$; the number of mixing component, $K$; size of
initialization random sample, $n_{0}$; the maximum number of iteration
$L_{0}$; the number of random starts $J$

Output: Partition $\mathcal{C^{*}}=\bigcup_{k=1}^{K}C_{k}$; robust FMGR
parameter estimate $\boldsymbol{\theta}^{*}=\boldsymbol{\theta}^{J_{0}}$;
outlier set $U^{*}=U^{J_{0}}$

1

2for  _$j=0,...,J$ _ do

3 Initialization: $U^{(old)}=\\{1,...,N\\};U^{(cur)}=\emptyset;L=0$

4 for  _k=1,…,K_ do

5 Draw a random sample of size $n_{0}$ from set $\\{1,...,N\\}$, indexed by
$I_{k}$

6 Run robust linear regression:
$(M_{k},\beta_{k},\sigma_{k})=:\textbf{RLM}(y_{I_{K}}\sim X_{I_{k}})$

7 Initialize posterior probability:
$p_{i}k=\mathcal{N}(y_{i}-x_{i}^{T}\boldsymbol{\beta}_{k};0,\sigma_{k}^{2}),i=1,...,N$

8

9 end for

10 while  _$U^{(old)}\neq U^{(cur)}\ &L<L_{0}$_ do

11 Let $U^{(old)}=U^{(cur)};L=L+1$

12 for  _k=1,…,K_ do

13 Let $I_{k}$ be sample indices most likely in cluster $k$

14 Let $U_{k}^{(cur)}$ be type 1 outliers of $Y_{I_{k}}\sim X_{I_{k},:}$,
using least trimmed suares robust regression

15

16 end for

17 $U^{(cur)}_{reg}=\bigcup_{k}U_{k}^{(cur)};S=\\{1,...,N\\}-U^{(cur)}_{reg}$

18 Update $(\boldsymbol{\theta},w)$ by HMR with the rest of samples in $S$,
let $U^{(cur)}_{spa}$ be type 2 outliers, and $W$ be the hybrid posterior
probability

19

20 end while

21
$U^{j}=U^{(cur)}_{reg}+U^{(cur)}_{spa},\boldsymbol{\theta}^{j}=\boldsymbol{\theta},w^{j}=w$

22 Let $F^{j}$ be a length-N binary vector whose $i-$th entry is 1 only if
$i\in U^{j}$

23 end for

Denote $J_{0}$ as the one such that $F_{J_{0}}$ is closet to the mean of
$\\{F_{j},j=1,...,J\\}$

Algorithm 2 Spatial Robust Mixture Regression (SRMR)

##  IV Experiments on synthetic data

We evaluated the performance of SRMR and selected baseline methods on a
comprehensive setup of synthetic datasets, and evaluated the overall
performance of SRMR in solving the spatially dependent mixture regression
problem with different number of mixture models, level of linear dependency,
spatial distribution, and ratio of outliers.

###  IV-A Baseline Methods

We collected in total nine existing methods to represent the current works. In
the field of mixture regression, Pan et al. [28] proposed DC-ADMM which
cluster mixture content in a group pursuit way. It has an implementation as
“PRclust” 111https://github.com/ChongWu-Biostat/prclust R package. In the
field of robust mixture regression, we collected two state-of-the-art
algorithms, Trimmed Likelihood Estimation (TLE) and Component-wise adaptive
Trimming Likelihood Estimation (CTLE) from R package “RobMixReg”
222https://cran.r-project.org/web/packages/RobMixReg/ in CRAN [27]. In the
field of spatial smooth regression, we collected four algorithms, spatially
clustered coefficient regression (SCC)333https://github.com/furong-
tamu/Supplementary-files-for-SCC [23], Spatialculster
444https://github.com/mpadge/spatialcluster [29], Spdep
555https://github.com/r-spatial/spdep/ [30], and ClustGeo
666https://cran.r-project.org/web/packages/ClustGeo/ [31]. However, only
ClustGeo can be executed under our formulation. In the field of segmentation
methods based on Markov Random Field, we collected two methods FRGMM
777https://sites.google.com/site/nguyen1j/home/10-code [26] and mrf2d
888https://freguglia.github.io/mrf2d/ [32]. However, these two methods aim to
clustering image pixels, which requires natural spatial orders from
neighborhood pixels as input, and hence cannot be applied to solve our
problem. Finally, we used four baseline methods DC-ADMM, TLE, CTLE, and
ClustGeo to perform comparison experiments. All baseline methods used with
their default parameters, except nit parameter in TLE and CTLERob were set as
10. For DC-ADMM, we used stability-prclust function to select the best
parameter, followed by the instruction. For ClustGeo, we used choicealpha
function to select best parameter.

###  IV-B Experimental setup

To simulate spatially dependent linear relationships, we first generate a
univariate independent variable $x$ from uniform distribution $X\sim U(-2,2)$
and dependent variable $y$ by
$y_{i}=\beta_{k}x_{i}+\sigma_{k},k=1,...,K,i=1,...,n$, where $K$ is number of
mixture models and $\beta$ is regression coefficient. Spatial coordinate of
each sample $s_{i}$ was generated from a multivariate normal distribution
$\mathcal{N}(\mu,\Sigma)$, where $\mu$ determines the center and $\Sigma$
determines the range and shape of each spatial region. We use
$\mu_{1}=[1,1]^{T},\mu_{2}=[-1,-1]^{T},\Sigma=diag(0.1,0.1)$ as the default
experimental setting, i.e. $K=2$ of two distinct and non-uniformly distributed
spatial regions. The two types of outliers wer further simulated. We simulated
Type 1 outliers by a rejection sampling approach. Specifically, we first
samples independent $(x,y)$ from $(U(-2,2),U(-8,8))$ and only accept the ones
whose Euclidean distance to the regression lines larger than two as Type 1
outliers. To simulate the Type 2 outliers of a certain ratio, we randomly
select the ratio of samples and reverse their spatial coordinate
$s_{i}=(c^{1}_{i},c^{2}_{i})$ by $s_{i}^{o}=(-c^{1}_{i},-c^{2}_{i})$.

We conducted the synthetic data based experiments for three types of method
evaluation:
(1) We evaluated the general performance of SRMR and baseline methods in
solving the spatially dependent mixture regression problem by the following
experimental setups (Fig 1A). Each time we perturbed one of the five factors
and fixed the others, including number of mixture regression models
$K=\\{2,3,4\\}$, total sample size $N=\\{100,200,400\\}$, error of linear
regression $\sigma=\\{0.1,0.2,0.5\\}$, rate of samples belong to (model1,
model2, outliers)=$\\{(0.4,0.4,0.2),(0.5,0.3,0.2),(0.6,0.2,0.2)\\}$ (only for
$K$=2), and coefficients of linear regression model
$\boldsymbol{\beta}=\\{(1.5,1.0),(1.5,0.1),(1.5,-1.2)\\}$ (only for $K$=2).
(2) We validated the robustness of SRMR and baseline methods in handling the
two types of outliers, namely Type 1 and 2 outliers by perturbing their ratio
from 10$\%$ to 20$\%$ (Fig 1B).
(3) We validated the capability of SRMR and baseline methods in detecting
different shapes and distributions of spatial regions. We simulated the
spatial coordinates from a multivariate normal distribution or a multivariate
uniform distribution, the former one simulates a round and dense spatial
region while the later one generates uniformly distributed 2D coordinates. The
simulated shapes are showcased in Fig 1C. In addition, we also evaluated if
SRMR is sensitive to different relative positions of the spatial regions. We
simulated two types of relative location of spatial regions, namely (i)
diagonal distribution by setting $\mu_{1}=[1,1]^{T},\mu_{2}=[-1,-1]^{T}$ and
(ii) horizontal distribution by setting
$\mu_{1}=[0.5,0]^{T},\mu_{2}=[-0.5,0]^{T}$. To simulate spatial regions of
imbalanced densities, we perturbed the covariance matrix of the spatial
coordinates from $diag(0.1,0.1)$ to $diag(0.5,0.1)$.

In summary, we set ten perturbation scenarios (Fig 1), each contains 2-3
different parameter settings. We conducted 100 replicates for each parameter
set in each scenario. In total, we obtained 2,500 synthetic data sets. The
mean value of evaluation metrics were used for performance evaluation.

Figure 1: Experiment Setting. Sub-figures without grid represent linear
relationship and sub-figures with grid represent spatial coordinates. For (b)
and (c), we only show partial plot which control factor is changed instead of
full plot (linear relationship and spatial coordinate) as (a). (a) contains
five different scenarios in terms of mixture regression. (b) contains two
scenarios to deal with Type 1 and Type 2 outliers. (c) contains three
scenarios for detecting different shapes and distributions of spatial regions.

###  IV-C Evaluation Metrics

We evaluated the performance of SRMR and baseline methods on synthetic
datasets, based on how accurate the methods can identify the simulated mixture
regression models and corresponding spatial regions, and distinguish the two
types of outliers. Four evaluation metrics were utilized in the synthetic data
based evaluations:

1) Rand Index (RI): $=\frac{\text{number of agreeing pairs}}{\text{number of
total pairs}}$ computes a similarity measure between two clusters by
considering counting the sample pairs that are assigned in the same or
different clusters in the predicted and true clusters.

2) Adjust Rand Index (ARI): $=\frac{\text{RI- Expected (RI)}}{\text{max(RI) -
Expected(RI)}}$, which is a corrected-for-chance version of RI.

3) Accuracy Rate (ACC) for outliear detection. ACC: $=\frac{\text{detected
true outliers}}{\text{true outliers}}$ measures the accuracy for
distinguishing the Type 1 and Type 2 outliers.

4) Error of Predicted Coefficients (PCE):
$=\sum_{k=1}^{K}(\beta_{k}-\beta^{p}_{l(k)})^{2}$ measures the distance
between the true regression coefficient $\beta_{k}$ of the regression
components $k=1,...,K$ and predicted regression coefficient $\beta^{p}$. Here
$l(k)=\underset{j}{\operatorname{argmin}}(\beta_{k}-\beta^{p}_{j})^{2}$, i.e.,
$\beta^{p}_{l}(k)$ is the predicted coefficient closest to $\beta_{k}$.

###  IV-D Performance

We organized the synthetic data experiment results in Table 1 into three
sections: mixture regression, robustness and spatial patterns. Overall, SRMR
outperforms baseline methods in all 10 experiment settings under almost all
evaluation metrics.

In Table 1, the first section (1st- 5th blocks) illustrated the performance of
SRMR and other methods in terms of the accuracies in detecting the
heterogeneous linear dependencies in different scenarios, with regards to
sample size, number of components, noise level, cluster balance and strength
of regression coefficients. SRMR could detect the clusters and regression
coefficients for each cluster very accurately, for different sample sizes,
components, and it is robust to the different noise levels, imbalance of
cluster sizes and small regression coefficients. Notably, because it
incorporates spatial information, it is able to differentiate two clusters
with very similar regression coefficients but different spatial locations.
Since DC-ADMM and ClustGeo are designed for clustering, but not regression,
the evaluation metrics ACC and PCE for these two methods are filled with NaN.
Although DC-ADMM proposed using a novel formulation for clustering, it cannot
handle outliers or incorporate spatial information. Thus, the performance of
DC-ADMM is the lowest in most of cases. As noise level of regression line
increased, the power of ordinary robust mixture regression methods TLE and
CTLE decreased, leading to lower RI and ARI score. When the clusters become
more and more imbalanced, the RI and ARI scores of all of the baseline methods
get much worse. When two clusters have very similar regression parameters, but
are far away in terms of spatial locations, TLE and CTLE cannot differentiate
the two clusters, as they didn’t account for spatial proximity, causing low RI
and ARI score.

The second section (6th-7th blocks) of Table 1 illustrated the performance of
all methods in terms of robustness to outlier contamination, including Type 1
outliers and Type 2 outliers. SRMR is highly robust to both regression
outliers and spatial outliers, and the clustering accuracies and parameter
estimates are almost unaffected in the presence of outliers. This is because
SRMR adopted a trimmed likelihood approach, and it is expected that the
outliers will not be taken into model estimations. Since DC-ADMM and ClustGeo
are not designed to handle the neither Type 1 or Type 2 outliers, their
performance consistently worse than TLE, CTLE, and SRMR. While TLE and CTLE
could handle regression outliers, they have no control over the spatial
proximity, and hence they are very sensitive Type 2 outliers, i.e., spatial
outliers. ACC of TLE and CTLE is around 70% due to spatial heterogeneity while
SRMR has 100% accuracy rate in all scenarios.

The third section (8th-10blocks) in Table 1 illustrate the performance of all
methods for different spatial patterns, regarding the shape, center and
density of the spatial clusters. SRMR is designed to detect heterogeneious
linear dependencies that is robust to both regression outliers and spatial
outliers, and its performance is consistently desirably with regards to
different spatial patterns. When the spatial distribution of the clusters are
changed from multivariate normal to multivariate uniform, it means the shape
of the clusters are less sphear, and more diffused. When the center of spatial
coordinate changed from diagonal to horizontal, the boundary of two spatial
centers became blurred, meaning there are more overlap between neighbouring
clusters. The performance of TLE, CTLE and ClustGeo got worse with more
cluster overlaps, while SRMR is robust to this complex situation thanks to the
integration of both regression and spatial similarity. ClustGEO is sensitive
to the imbalanced density of different clusters, while SRMR is unaffected.

In summary, SRMR is the only method that could model the linear dependency
between response and predictors that vary in the spatial domain, and detect
clusters of observations with both similarities in regression parameters and
spatial proximity. And it is robust to both outliers in regression fitting and
spatial locations. It has produced highly favorable performance in different
simulation settings, with regards to different levels of regression/spatial
noise, outliers, and mixture imbalance.

##  V Experiments on real-world data

We further validated SRMR on two real-world datasets, namely (1) a geospatial
economics data collected from 298 cities of China and (2) a spatial
transcriptomics data collected from 3,798 spatial spots on a 2D breast cancer
tissue. The synthetic data based experiments clearly suggested that SRMR is
the only method that can effectively solve the spatially dependent mixture
regression problem compared to the baseline methods. In the real-world data
based experiments, we mainly focused on illustrating the contextual meaning of
the spatial regions and corresponding regression models identified by SRMR. We
also evaluated the goodness of fitting and significance of the spatially
dependent mixture regression models as well as the running time of the tested
methods.

Figure 2: Real-world data based experiments. a1: SRMR, a2: SRMR, a3: TLE, a4:
CTLE, a5: DC-ADMM, a6: ClustGeo; a1,a3-a6: $ES\sim GDP$, a2: $Income\sim GDP$.
a1-a6: cities of different regression components are red, blue, green and
orange colored, while the outliers are colored by grey.

###  V-A Application on Geospatial Economics Data

We collected 7 economic features, namely total GDP, public income, public
spend, educational spend, technology spend, population, and averaged personal
income, and latitude and longitudes, for 298 cities in China. We evaluated
SRMR and baseline methods to this data set. We utilized each of the eight
features as a dependent variable and selected others as independent variables
When applying SRMR and other regression models, while all the features were
utilized as the input of ClustGeo. Similar to the synthetic data based
experiments, SRMR is the only method can identify spatially dependent mixture
regression models. In contrast, TLE, CTLERob, and DC-ADMM only detected
spatial independent regression models, and ClustGeo output a spatial
segmentation based on all features.

For a clear visualization and explanation, we illustrated two univariate
regressions of $Educational\ Spend\ (ES)\ \sim\ GDP$ and $Income\ \sim\ GDP$.
For both $ES$ and $Income$, SRMR identified four spatial regions corresponding
to the north-east, middle-east, south-east and west regions of China (Fig
2a1,2a2). The spatial regions detected by SRMR show distinct different
dependency of $ES$ and $Income$ with $GDP$. Specifically, $ES$ is positively
associated with $GDP$ in the middle-east ($ES=0.24\cdot GDP+10.9$) and north-
east China ($ES=0.4\cdot GDP+9.17$). The south-east cities have more stable
$ES$, which less depends on $GDP$ ($ES=0.8\cdot GDP+4.19$), while a negative
association of $ES$ and $GDP$ are observed in the west cities ($ES=-0.39\cdot
GDP+17.09$). The high dependency in middle-east and north-east cities and less
dependency in south-east cities are consistent to our knowledge, as the
middle-east and north-east China are promoting the education system basis
while the education systems south-east China are relatively stable. We also
checked the cities in the west China that have high $GDP$ but low $ES$. Such
cities include Dongying, Ordos, Karamay, etc., which are developing more neo
energy business rather than education in the recent years. Similar
observations were also made in the $Income\ \sim\ GDP$ model (Fig 2a2). The
SRMR outputs suggested the personal $Income$ in the north-east, south-east and
west cities less depends on $GDP$ while more positive dependency between
$Income$ and $GDP$ was observed in middle-east cities, especially the well
developed cities Beijing, Shanghai, Tianjin, Hangzhou, etc. On the other hand,
on both $Educational\ Spend\ (ES)\ \sim\ GDP$ and $Income\ \sim\ GDP$, TLE and
CTLERob failed to identify such spatial dependent and contextual meaningful
patterns while both of them tend to over-fit the mixture of regressions (Fig
2a3, 2a4). DC-ADMM identified all cities as one class (Fig 2a5) while ClustGeo
identified three distinct non-overlapping spatial regions without offering
explainable regional specific feature dependencies (Fig 2a6).

###  V-B Application on Spatial Transcriptomics Data

10x Genomics spatial transcriptomics (ST) is a recent commercialized technique
to measure spatial coordinates associated gene expression signal from a
biological tissue sample, which it has a broad utilization in biomedical
research. A typical ST data is a matrix consisting of $\sim$15,000 genes
(rows) in $\sim$4,000 individual spatial spots (columns), and each spot has a
2D spatial coordinate (Fig 2b1). The spatial spots are uniformly distributed.
A key challenge in ST data analysis is to infer the spatially dependent and
biologically meaningful functional variations from the high dimensional
feature matrix (genes by spatial spots). Here we illustrate that SRMR enables
a new type of ST data analysis by simultaneously identify spatial regions in
which the expression level of genes show different level of dependency, which
directly annotate the biological meaning of each detected region.

We applied SRMR and baseline methods on the v1.1 ST data of breast cancer
provided by 10xgenomics.com, consisting of 13,161 genes and 3,798 spatial
spots. We first selected 500 genes that having high expression level and
having known tumor micro-environment related functions. We fit the regression
model Gene${}_{1}\ \sim$ Gene2 for each pair of the 500 genes by using SRMR,
TLE, CTLERob and DC-ADMM and conducted ClustGeo by using all the 500 genes.
Similar to the synthetic and Geospatial data, SRMR is the only method that
detected spatially dependent mixture regression models in the ST data. General
spatial segmentation, such as ClustGeo, identifies spatial regions by using
the whole feature matrix (Fig 2b3), which is consistent to the distribution of
the averaged gene expression signal level (Fig 2b2). On the other hand, we
identified more than 500 overlapped spatial regions by using SRMR, each having
varied dependency among certain genes. Fig 2b4 showcased two distinct spatial
regions only identified by SRMR, which have varied dependency between the
CD79A and CD79B genes as shown in Fig 2b5. CD79A/B are key genes involved in
maturation and functional variation of B cells. The varied dependency of CD79A
and CD79B characterizes distinct sub-regions in one breast cancer tissue that
potentially have different immune activities and responses to immuno-therapy.

In summary, compared with baseline methods, SRMR is the only method can
effectively solve the spatially dependent mixture regression problem on the
two real-world data. For the analysis of a single regression model in the
real-world data, the running time of SRMR, TLE and CTLERob are about 15s, 10s
and 2s, respectively. The running time of SRMR is slower, but also comparable
to the baseline robust mixture regression approaches. The running time of DC-
ADMM and ClustGeo are about 0.01s.

##  VI Conclusion

We developed a new statistical model of high dimensional data with matched
spatial information, namely spatially dependent mixture regression. We also
developed spatial robust mixture regression (SRMR) analysis as an effective
solution of the problem. SRMR is empowered by an inference scheme to assess
statistical significance of spatial dependent finite mixture regression
models. On both synthetic and real-world data based experiments, we
demonstrated that SRMR is the only capability can solve the spatially
dependent mixture regression problem. Particularly, SRMR enables a new type of
spatial segmentation analysis by detecting large sets of spatial regions
having varied dependency among certain features. Compared with conventional
spatial segmentation analysis, the regions identified by SRMR characterize
more spatial dependent variations conceived in the data and enable better
contextual explanation. The source codes of SRMR and the analysis of this
study are provided at https://github.com/changwn/SRMR.

## References

  * Dolde and Tirtiroglu [1997] W. Dolde and D. Tirtiroglu, “Temporal and spatial information diffusion in real estate price changes and variances,” _Real Estate Economics_ , vol. 25, no. 4, pp. 539–565, 1997.
  * Van Den Heuvel and Pol [2010] M. P. Van Den Heuvel and H. E. H. Pol, “Exploring the brain network: a review on resting-state fmri functional connectivity,” _European neuropsychopharmacology_ , vol. 20, no. 8, pp. 519–534, 2010.
  * McLachlan and Peel [2004] G. J. McLachlan and D. Peel, _Finite mixture models_. John Wiley & Sons, 2004.
  * Frühwirth-Schnatter [2006] S. Frühwirth-Schnatter, _Finite mixture and Markov switching models_. Springer Science & Business Media, 2006.
  * Yu et al. [2020] C. Yu, W. Yao, and G. Yang, “A selective overview and comparison of robust mixture regression estimators,” _International Statistical Review_ , vol. 88, no. 1, pp. 176–202, 2020.
  * Markatou [2000] M. Markatou, “Mixture models, robustness, and the weighted likelihood methodology,” _Biometrics_ , vol. 56, no. 2, pp. 483–486, 2000.
  * Shen et al. [2004] H.-b. Shen, J. Yang, and S.-t. Wang, “Outlier detecting in fuzzy switching regression models,” in _International Conference on Artificial Intelligence: Methodology, Systems, and Applications_. Springer, 2004, pp. 208–215.
  * Bai et al. [2012] X. Bai, W. Yao, and J. E. Boyer, “Robust fitting of mixture regression models,” _Computational Statistics & Data Analysis_, vol. 56, no. 7, pp. 2347–2359, 2012.
  * Bashir and Carter [2012] S. Bashir and E. Carter, “Robust mixture of linear regression models,” _Communications in Statistics-Theory and Methods_ , vol. 41, no. 18, pp. 3371–3388, 2012.
  * Song et al. [2014] W. Song, W. Yao, and Y. Xing, “Robust mixture regression model fitting by laplace distribution,” _Computational Statistics & Data Analysis_, vol. 71, pp. 128–137, 2014.
  * Yao et al. [2014] W. Yao, Y. Wei, and C. Yu, “Robust mixture regression using the t-distribution,” _Computational Statistics & Data Analysis_, vol. 71, pp. 116–127, 2014.
  * Peel and McLachlan [2000] D. Peel and G. J. McLachlan, “Robust mixture modelling using the t distribution,” _Statistics and computing_ , vol. 10, no. 4, pp. 339–348, 2000.
  * Yu et al. [2017] C. Yu, W. Yao, and K. Chen, “A new method for robust mixture regression,” _Canadian Journal of Statistics_ , vol. 45, no. 1, pp. 77–94, 2017.
  * Neykov et al. [2007] N. Neykov, P. Filzmoser, R. Dimova, and P. Neytchev, “Robust fitting of mixtures using the trimmed likelihood estimator,” _Computational Statistics & Data Analysis_, vol. 52, no. 1, pp. 299–308, 2007.
  * Doğru and Arslan [2018] F. Z. Doğru and O. Arslan, “Robust mixture regression modeling using the least trimmed squares (lts)-estimation method,” _Communications in Statistics-Simulation and Computation_ , vol. 47, no. 7, pp. 2184–2196, 2018.
  * García-Escudero et al. [2010] L. A. García-Escudero, A. Gordaliza, A. Mayo-Íscar, and R. San Martín, “Robust clusterwise linear regression through trimming,” _Computational Statistics & Data Analysis_, vol. 54, no. 12, pp. 3057–3069, 2010.
  * García-Escudero et al. [2017] L. A. García-Escudero, A. Gordaliza, F. Greselin, S. Ingrassia, and A. Mayo-Íscar, “Robust estimation of mixtures of regressions with random covariates, via trimming and constraints,” _Statistics and Computing_ , vol. 27, no. 2, pp. 377–402, 2017.
  * Brunsdon et al. [1996] C. Brunsdon, A. S. Fotheringham, and M. E. Charlton, “Geographically weighted regression: a method for exploring spatial nonstationarity,” _Geographical analysis_ , vol. 28, no. 4, pp. 281–298, 1996.
  * Fotheringham et al. [2003] A. S. Fotheringham, C. Brunsdon, and M. Charlton, _Geographically weighted regression: the analysis of spatially varying relationships_. John Wiley & Sons, 2003.
  * Wheeler and Páez [2010] D. C. Wheeler and A. Páez, “Geographically weighted regression,” in _Handbook of applied spatial analysis_. Springer, 2010, pp. 461–486.
  * Fuentes [2002] M. Fuentes, “Spectral methods for nonstationary spatial processes,” _Biometrika_ , vol. 89, no. 1, pp. 197–210, 2002.
  * Banerjee et al. [2014] S. Banerjee, B. P. Carlin, and A. E. Gelfand, _Hierarchical modeling and analysis for spatial data_. CRC press, 2014\.
  * Li and Sang [2019] F. Li and H. Sang, “Spatial homogeneity pursuit of regression coefficients for large datasets,” _Journal of the American Statistical Association_ , vol. 114, no. 527, pp. 1050–1062, 2019.
  * Tibshirani et al. [2005] R. Tibshirani, M. Saunders, S. Rosset, J. Zhu, and K. Knight, “Sparsity and smoothness via the fused lasso,” _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , vol. 67, no. 1, pp. 91–108, 2005\.
  * Nguyen and Wu [2011] T. M. Nguyen and Q. J. Wu, “Gaussian-mixture-model-based spatial neighborhood relationships for pixel labeling problem,” _IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics)_ , vol. 42, no. 1, pp. 193–202, 2011.
  * Nguyen and Wu [2012] ——, “Fast and robust spatially constrained gaussian mixture model for image segmentation,” _IEEE transactions on circuits and systems for video technology_ , vol. 23, no. 4, pp. 621–635, 2012.
  * Chang et al. [2020] W. Chang, C. Wan, C. Yu, W. Yao, C. Zhang, and S. Cao, “Robmixreg: an r package for robust, flexible and high dimensional mixture regression,” _bioRxiv_ , 2020.
  * Wu et al. [2016] C. Wu, S. Kwon, X. Shen, and W. Pan, “A new algorithm and theory for penalized regression-based clustering,” _The Journal of Machine Learning Research_ , vol. 17, no. 1, pp. 6479–6503, 2016.
  * Guo [2008] D. Guo, “Regionalization with dynamically constrained agglomerative clustering and partitioning (redcap),” _International Journal of Geographical Information Science_ , vol. 22, no. 7, pp. 801–823, 2008.
  * Bivand et al. [2011] R. Bivand, L. Anselin, O. Berke, A. Bernat, M. Carvalho, Y. Chun, C. Dormann, S. Dray, R. Halbersma, N. Lewin-Koh _et al._ , “spdep: Spatial dependence: weighting schemes, statistics and models,” 2011.
  * Chavent et al. [2018] M. Chavent, V. Kuentz-Simonet, A. Labenne, and J. Saracco, “Clustgeo: an r package for hierarchical clustering with spatial constraints,” _Computational Statistics_ , vol. 33, no. 4, pp. 1799–1822, 2018.
  * Freguglia and Garcia [2020] V. Freguglia and N. L. Garcia, “mrf2d: Markov random field image models in r,” _arXiv preprint arXiv:2006.00383_ , 2020.


11institutetext: 1 FIT BUT, Božetěchova 2, 61266 Brno, Czech Republic
2 VERIMAG, UJF, 2. av. de Vignate, 38610 Gières, France
email: 11email<EMAIL_ADDRESS>
# Optimizing an LTS-Simulation Algorithm††thanks: This is a version of
Technical Report No. FIT-TR-2009-03 of Faculty of Information Technology, Brno
University of Technology, accompanying the work published originally as [10].

Lukáš Holík1 Jiří Šimáček1,2

###### Abstract

When comparing the fastest algorithm for computing the largest simulation
preorder over Kripke structures with the one for labeled transition systems
(LTS), there is a noticeable time and space complexity blow-up proportional to
the size of the alphabet of an LTS. In this paper, we present optimizations
that suppress this increase of complexity and may turn a large alphabet of an
LTS to an advantage. Our experimental results show significant speed-ups and
memory savings. Moreover, the optimized algorithm allows one to improve
asymptotic complexity of procedures for computing simulations over tree
automata using recently proposed algorithms based on computing simulation over
certain special LTS derived from a tree automaton.

##  1 Introduction

A practical limitation of automated methods dealing with LTSs—such as LTL
model checking, regular model checking, etc.—is often the size of generated
LTSs. One of the well established approaches to overcome this problem is the
reduction of an LTS using a suitable equivalence relation according to which
the states of the LTS are collapsed. A good candidate for such a relation is
simulation equivalence. It strongly preserves logics like $ACTL^{*}$,
$ECTL^{*}$, and $LTL$ [6, 7, 8], and with respect to its reduction power and
computation cost, it offers a desirable compromise among the other common
candidates, such as bisimulation equivalence [11, 13] and language
equivalence. The currently fastest LTS-simulation algorithm (denoted as
LRT—labeled RT) has been published in [1]. It is a straightforward
modification of the fastest algorithm (denoted as RT, standing for Ranzato-
Tapparo) for computing simulation over Kripke structures [12], which improves
the algorithm from [8]. The time complexity of RT amounts to
$\mathcal{O}(|P_{\mathit{Sim}}||\delta|)$, space complexity to
$\mathcal{O}(|P_{\mathit{Sim}}||S|)$. In the case of LRT, we obtain time
complexity
$\mathcal{O}(|P_{\mathit{Sim}}||\delta|+|\Sigma||P_{\mathit{Sim}}||S|)$ and
space complexity $\mathcal{O}(|\Sigma||P_{\mathit{Sim}}||S|)$. Here, $S$ is
the set of states, $\delta$ is the transition relation, $\Sigma$ is the
alphabet and $P_{\mathit{Sim}}$ is the partition of $S$ according to the
simulation equivalence. The space complexity blow-up of LRT is caused by
indexing the data structures of RT by the symbols of the alphabet.

In this paper, we propose an optimized version of LRT (denoted OLRT) that
lowers the above described blow-up. We exploit the fact that not all states of
an LTS have incoming and outgoing transitions labeled by all symbols of an
alphabet, which allows us to reduce the memory footprint of the data
structures used during the computation. Our experiments show that the
optimizations we propose lead to significant savings of space as well as of
time in many practical cases. Moreover, we have achieved a promising reduction
of the asymptotic complexity of algorithms for computing tree-automata
simulations from [1] using OLRT.

##  2 Preliminaries

Given a binary relation $\rho$ over a set $X$, we use $\rho(x)$ to denote the
set $\\{y\mid(x,y)\in\rho\\}$. Then, for a set $Y\subseteq X$,
$\rho(Y)=\bigcup\\{\rho(y)\mid y\in Y\\}$. A _partition-relation pair_ over
$X$ is a pair $\langle P,Rel\rangle$ where $P\subseteq 2^{X}$ is a partition
of $X$ (we call elements of $P$ _blocks_) and $Rel\subseteq P\times P$. A
partition-relation pair $\langle P,Rel\rangle$ _induces_ the relation
${\rho}={\bigcup_{(B,C)\in Rel}B\times C}$. We say that $\langle
P,\mathit{Rel}\rangle$ is the _coarsest partition-relation pair_ inducing
$\rho$ if any two $x,y\in X$ are in the same block of $P$ if and only if
$\rho(x)=\rho(y)$ and $\rho^{-1}(x)=\rho^{-1}(y)$. Note that in the case when
$\rho$ is a preorder and $\langle P,\mathit{Rel}\rangle$ is coarsest, then $P$
is the set of equivalence classes of $\rho\cap\rho^{-1}$ and $\mathit{Rel}$ is
a partial order.

A _labeled transition system (LTS)_ is a tuple $T=(S,\Sigma,\\{\delta_{a}\mid
a\in\Sigma\\})$, where $S$ is a finite set of states, $\Sigma$ is a finite set
of labels, and for each $a\in\Sigma$, $\delta_{a}\subseteq{S\times S}$ is an
$a$-labeled transition relation. We use $\delta$ to denote
$\bigcup_{a\in\Sigma}\delta_{a}$. A simulation over $T$ is a binary relation
$\rho$ on $S$ such that if $(u,v)\in\rho$, then for all $a\in\Sigma$ and
$u^{\prime}\in\mathrm{\delta}_{a}(u)$, there exists
$v^{\prime}\in\mathrm{\delta}_{a}(v)$ such that
$(u^{\prime},v^{\prime})\in\rho$. It can be shown that for a given LTS $T$ and
an initial preorder $\mathit{I}\subseteq S\times S$, there is a unique maximal
simulation ${\mathit{Sim}_{\mathit{I}}}$ on $T$ that is a subset of
$\mathit{I}$, and that ${\mathit{Sim}_{\mathit{I}}}$ is a preorder (see [1]).

##  3 The Original LRT Algorithm

In this section, we describe the original version of the algorithm presented
in [1], which we denote as LRT (see Algortihm 1).

Input: an LTS $T=(S,\Sigma,\\{\delta_{a}\mid a\in\Sigma\\})$, partition-
relation pair $\langle P_{\mathit{I}},\mathit{Rel}_{\mathit{I}}\rangle$

Output: partition-relation pair $\langle P,\mathit{Rel}\rangle$

1

/* initialization */

2 $\langle P,\mathit{Rel}\rangle\leftarrow\langle
P_{I},\mathit{Rel}_{I}\rangle$ /* $\leftarrow\langle
P_{\mathit{I}\cap\mathit{Out}},\mathit{Rel}_{\mathit{I}\cap\mathit{Out}}\rangle$
*/

3forall  _$B\in P$ and $a\in\Sigma$_ do /* $a\in\mathrm{in}(B)$ */

4 forall  _$v\in S$_ do
$\mathit{Count}_{a}(v,B)=|\mathrm{\delta}_{a}(v)\cap\bigcup\mathit{Rel}(B)|$ ;

5 /* $v\in\mathrm{\delta}^{-1}_{a}(S)$ */ $\mathit{Remove}_{a}(B)\leftarrow
S\setminus\mathrm{\delta}^{-1}_{a}(\bigcup\mathit{Rel}(B))$ /*
$\leftarrow\mathrm{\delta}^{-1}_{a}(S){\setminus}\mathrm{\delta}^{-1}_{a}(\bigcup\mathit{Rel}(B))$
*/

/* computation */

6 while  _exists $B\in P$ and $a\in\Sigma$ such that
$\mathit{Remove}_{a}(B)\neq\emptyset$_ do

7 $\mathit{Remove}\leftarrow\mathit{Remove}_{a}(B)$;

8 $\mathit{Remove}_{a}(B)\leftarrow\emptyset$;

9 $\langle
P_{\mathsf{prev}},\mathit{Rel}_{\mathsf{prev}}\rangle\leftarrow\langle
P,\mathit{Rel}\rangle$;

10 $P\leftarrow\mathit{Split}(P,\mathit{Remove})$;

11 $\mathit{Rel}\leftarrow\\{(C,D)\in P\times
P\mid(C_{\mathsf{prev}},D_{\mathsf{prev}})\in\mathit{Rel}_{\mathsf{prev}}\\}$;

12 forall  _$C\in P$ and $b\in\Sigma$_ do /* $b\in\mathrm{in}(C)$ */

13 $\mathit{Remove}_{b}(C)\leftarrow\mathit{Remove}_{b}(C_{\mathsf{prev}})$;

14 forall  _$v\in S$_ do
$\mathit{Count}_{b}(v,C)\leftarrow\mathit{Count}_{b}(v,C_{\mathsf{prev}})$ ;

15 /* $v\in\mathrm{\delta}^{-1}_{b}(S)$ */

16 forall  _$C\in P$ such that
$C\cap\mathrm{\delta}^{-1}_{a}(B)\neq\emptyset$_ do

17 forall  _$D\in P$ such that $D\subseteq\mathit{Remove}$_ do

18 if  _$(C,D)\in\mathit{Rel}$_ then

19 $\mathit{Rel}\leftarrow\mathit{Rel}\setminus\\{(C,D)\\}$;

20 forall  _$b\in\Sigma$ and $v\in\mathrm{\delta}^{-1}_{b}(D)$_ do /*
$b\in\mathrm{in}(D)\cap\mathrm{in}(C)$ */

21 $\mathit{Count}_{b}(v,C)\leftarrow\mathit{Count}_{b}(v,C)-1$;

22 if  _$\mathit{Count}_{b}(v,C)=0$_ then
$\mathit{Remove}_{b}(C)\leftarrow\mathit{Remove}_{b}(C)\cup\\{v\\}$ ;

23

Algorithm 1 LRT Algorithm

The algorithm gradually refines a partition-relation pair $\langle
P,\mathit{Rel}\rangle$, which is initialized as the coarsest partition-
relation pair inducing an initial preorder $\mathit{I}$. After its
termination, $\langle P,\mathit{Rel}\rangle$ is the coarsest partition-
relation pair inducing ${\mathit{Sim}_{\mathit{I}}}$. The basic invariant of
the algorithm is that the relation induced by $\langle P,\mathit{Rel}\rangle$
is always a superset of ${\mathit{Sim}_{\mathit{I}}}$.

The while-loop refines the partition $P$ and then prunes the relation
$\mathit{Rel}$ in each iteration of the while-loop. The role of the
$\mathit{Remove}$ sets can be explained as follows: During the initialization,
every $\mathit{Remove}_{a}(B)$ is filled by states $v$ such that
$\mathrm{\delta}_{a}(v)\cap\bigcup\mathit{Rel}(B)=\emptyset$ (there is no
$a$-transition leading from $v$ “above” $B$ wrt. $\mathit{Rel}$). During the
computation phase, $v$ is added into $\mathit{Remove}_{a}(B)$ after
$\mathrm{\delta}_{a}(v)\cap\bigcup\mathit{Rel}(B)$ becomes empty (because of
pruning $\mathit{Rel}$ on line 17). Emptiness of
$\mathrm{\delta}_{a}(v)\cap\mathit{Rel}(B)$ is tested on line 20 using
counters $\mathit{Count}_{a}(v,B)$, which record the cardinality of
$\mathrm{\delta}_{a}(v)\cap\mathit{Rel}(B)$. From the definition of
simulation, and because the relation induced by $\langle
P,\mathit{Rel}\rangle$ is always a superset of ${\mathit{Sim}_{\mathit{I}}}$,
$\mathrm{\delta}_{a}(v)\cap\bigcup\mathit{Rel}(B)=\emptyset$ implies that for
all $u\in\mathrm{\delta}^{-1}_{a}(B)$,
$(u,v)\not\in{\mathit{Sim}_{\mathit{I}}}$ ($v$ cannot simulate any
$u\in\mathrm{\delta}^{-1}_{a}(B)$). To reflect this, the relation
$\mathit{Rel}$ is pruned each time $\mathit{Remove}_{a}(B)$ is processed. The
code on lines 8–13 prepares the partition-relation pair and all the data
structures. First, $\mathit{Split}(P,\mathit{Remove}_{a}(B))$ divides every
block $B^{\prime}$ into $B^{\prime}\cap\mathit{Remove}_{a}(B)$ (which cannot
simulate states from $\mathrm{\delta}^{-1}_{a}(B)$ as they have empty
intersection with $\mathrm{\delta}^{-1}_{a}(\mathit{Rel}(B))$), and
$B^{\prime}\setminus\mathit{Remove}_{a}(B)$. More specifically, for a set
$\mathit{Remove}\subseteq S$, $\mathit{Split}(P,\mathit{Remove})$ returns a
finer partition $P^{\prime}=\\{B\setminus\mathit{Remove}\mid B\in
P\\}\cup\\{B\cap\mathit{Remove}\mid B\in P\\}$. After refining $P$ by the
$\mathit{Split}$ operation, the newly created blocks of $P$ inherit the data
structures (counters $\mathit{Count}$ and $\mathit{Remove}$ sets) from their
“parents” (for a block $B\in P$, its parent is the block $B_{\mathsf{prev}}\in
P_{\mathsf{prev}}$ such that $B\subseteq B_{\mathsf{prev}}$). $\mathit{Rel}$
is then updated on line 17 by removing the pairs $(C,D)$ such that
$C\cap\mathrm{\delta}^{-1}_{a}(B)\neq\emptyset$ and
$D\subseteq\mathit{Remove}_{a}(B)$. The change of $\mathit{Rel}$ causes that
for some states $u\in S$ and symbols $b\in\Sigma$,
$\mathrm{\delta}_{a}(u)\cap\bigcup\mathit{Rel}(C)$ becomes empty. To propagate
the change of the relation along the transition relation, $u$ will be moved
into $\mathit{Remove}_{b}(C)$ on line 20, which will cause new changes of the
relation in the following iterations of the while-loop. If there is no
nonempty $\mathit{Remove}$ set, then $\langle P,\mathit{Rel}\rangle$ is the
coarsest partition-relation pair inducing ${\mathit{Sim}_{\mathit{I}}}$ and
the algorithm terminates. Correctness of LRT is stated by Theorem 3.1.

######  Theorem 3.1 ([1])

With an LTS $T=(S,\Sigma,\\{\delta_{a}\mid a\in\Sigma\\})$ and the coarsest
partition-relation pair $\langle
P_{\mathit{I}},\mathit{Rel}_{\mathit{I}}\rangle$ inducing a preorder
$I\subseteq S\times S$ on the input, LRT terminates with the coarsest
partition-relation pair $\langle P,\mathit{Rel}\rangle$ inducing
${\mathit{Sim}_{\mathit{I}}}$.

##  4 Optimizations of LRT

The optimization we are now going to propose reduces the number of counters
and the number and the size of $\mathit{Remove}$ sets. The changes required by
OLRT are indicated in Algorithm 1 on the right hand sides of the concerned
lines.

We will need the following notation. For a state $v\in S$,
$\mathrm{in}(v)=\\{a\in\Sigma\mid\mathrm{\delta}^{-1}_{a}(v)\neq\emptyset\\}$
is the set of _input symbols_ and
$\mathrm{out}(v)=\\{a\in\Sigma\mid\mathrm{\delta}_{a}(v)\neq\emptyset\\}$ is
the set of _output symbols_ of $v$. The _output preorder_ is the relation
$\mathit{Out}=\bigcap_{a\in\Sigma}\mathrm{\delta}^{-1}_{a}(S)\times\mathrm{\delta}^{-1}_{a}(S)$
(this is, ${(u,v)\in\mathit{Out}}$ if and only if
${\mathrm{out}(u)\subseteq\mathrm{out}(v)}$).

To make our optimization possible, we have to initialize $\langle
P,\mathit{Rel}\rangle$ by the finer partition-relation pair $\langle
P_{\mathit{I}\cap\mathit{Out}},\mathit{Rel}_{\mathit{I}\cap\mathit{Out}}\rangle$
(instead of $\langle P_{\mathit{I}},\mathit{Rel}_{\mathit{I}}\rangle$), which
is the coarsest partition-relation pair inducing the relation
$\mathit{I}\cap\mathit{Out}$. As both $\mathit{I}$ and $\mathit{Out}$ are
preorders, $\mathit{I}\cap\mathit{Out}$ is a preorder too. As
${\mathit{Sim}_{\mathit{I}}}\subseteq\mathit{I}$ and
${\mathit{Sim}_{\mathit{I}}}\subseteq\mathit{Out}$ (any simulation on $T$ is a
subset of $\mathit{Out}$), ${\mathit{Sim}_{\mathit{I}}}$ equals the maximal
simulation included in $\mathit{I}\cap\mathit{Out}$. Thus, this step itself
does not influence the output of the algorithm.

Assuming that $\langle P,\mathit{Rel}\rangle$ is initialized to $\langle
P_{\mathit{I}\cap\mathit{Out}},\mathit{Rel}_{\mathit{I}\cap\mathit{Out}}\rangle$,
we can observe that for any $B\in P$ and $a\in\Sigma$ chosen on line 5, the
following two claims hold:

###### Claim

If $a\not\in\mathrm{in}(B)$, then skipping this iteration of the while-loop
does not affect the output of the algorithm.

###### Proof.

In an iteration of the while-loop processing $\mathit{Remove}_{a}(B)$ with
$a\not\in\mathrm{in}(B)$, as there is no $C\in P$ with
$\delta_{a}(C)\cap\mathit{Rel}(B)\neq\emptyset$, the for-loop on line 16 stops
immediately. No pair $(C,D)$ will be removed from $\mathit{Rel}$ on line 17,
no counter will be decremented, and no state will be added into a
$\mathit{Remove}$ set. The only thing that can happen is that
$\mathit{Split}(P,\mathit{Remove})$ refines $P$. However, in this case, this
refinement of $P$ would be done anyway in other iterations of the while-loop
when processing sets $\mathit{Remove}_{b}(C)$ with $b\in\mathrm{in}(C)$. To
see this, note that correctness of the algorithm does not depend on the order
in which nonempty $\mathit{Remove}$ sets are processed. Therefore, we can
postpone processing all the nonempty $\mathit{Remove}_{a}(B)$ sets with
$a\not\in\mathrm{in}(B)$ to the end of the computation. Recall that processing
no of these $\mathit{Remove}$ sets can cause that an empty $\mathit{Remove}$
set becomes nonempty. Thus, the algorithm terminates after processing the last
of the postponed $\mathit{Remove}_{a}(B)$ sets. If processing some of these
$\mathit{Remove}_{a}(B)$ with $a\not\in\mathrm{in}(B)$ refines $P$, $P$ will
contain blocks $C,D$ such that both $(C,D)$ and $(D,C)$ are in $\mathit{Rel}$
(recall that when processing $\mathit{Remove}_{a}(B)$, no pair of blocks can
be removed from $\mathit{Rel}$ on line 17). This means that the final $\langle
P,\mathit{Rel}\rangle$ will not be coarsest, which is a contradiction with
Theorem 3.1. Thus, processing the postponed $\mathit{Remove}_{a}(B)$ sets can
influence nor $\mathit{Rel}$ neither $P$, and therefore they do not have to be
processed at all. ∎

###### Claim

It does not matter whether we assign $\mathit{Remove}_{a}(B)$ or
$\mathit{Remove}_{a}(B)\setminus(S\setminus\mathrm{\delta}^{-1}_{a}(S))$ to
$\mathit{Remove}$ on line 6.

###### Proof.

Observe that $v$ with $a\not\in\mathrm{out}(v)$ (i.e., $v\in
S\setminus\mathrm{\delta}^{-1}_{a}(S)$) cannot be added into
$\mathit{Remove}_{a}(B)$ on line 20, as this would mean that $v$ has an
$a$-transition leading to $D$. Therefore, $v$ can get into
$\mathit{Remove}_{a}(B)$ only during initialization on line 4 together with
all states from $S\setminus\mathrm{\delta}^{-1}_{a}(S)$. After
$\mathit{Remove}_{a}(B)$ is processed (and emptied) for the first time, no
state from $S\setminus\mathrm{\delta}^{-1}_{a}(S)$ can appear there again.
Thus, $\mathit{Remove}_{a}(B)$ contains states from
$S\setminus\mathrm{\delta}^{-1}_{a}(S)$ only when it is processed for the
first time and then it contains all of them. It can be shown that for any
partition $Q$ of a set $X$ and any $Y\subseteq X$, if $\mathit{Split}(Q,Y)=Q$,
then also for any $Z\subseteq X$ with $Y\subseteq Z$,
$\mathit{Split}(Q,Z)=\mathit{Split}(Q,Z\setminus Y)$. As $P$ refines
$P_{\mathit{I}\cap\mathit{Out}}$,
$\mathit{Split}(P,S\setminus\mathrm{\delta}^{-1}_{a}(S))=P$. Therefore, as
$S\setminus\mathrm{\delta}^{-1}_{a}(S)\subseteq\mathit{Remove}_{a}(B)$,
$\mathit{Split}(P,\mathit{Remove}_{a}(B))=\mathit{Split}(P,\mathit{Remove}_{a}(B)\setminus(S\setminus\mathrm{\delta}^{-1}_{a}(S)))$.
We have shown that removing $S\setminus\mathrm{\delta}^{-1}_{a}(S)$ from
$\mathit{Remove}$ does not influence the result of the $\mathit{Split}$
operation in this iteration of the while-loop (note that this implies that all
blocks from the new partition are included in or have empty intersection with
$S\setminus\mathrm{\delta}^{-1}_{a}(S)$). It remains to show that it also does
not influence updating $\mathit{Rel}$ on line 17. Removing
$S\setminus\mathrm{\delta}^{-1}_{a}(S)$ from $\mathit{Remove}$ could only
cause that the blocks $D$ such that $D\subseteq
S\setminus\mathrm{\delta}^{-1}_{a}(S)$ that were chosen on line 15 with the
original value of $\mathit{Remove}$ will not be chosen with the restricted
$\mathit{Remove}$. Thus, some of the pairs $(C,D)$ removed from $\mathit{Rel}$
with the original version of $\mathit{Remove}$ could stay in $\mathit{Rel}$
with the restricted version of $\mathit{Remove}$. However, such a pair $(C,D)$
cannot exist because with the original value of $\mathit{Remove}$, if $(C,D)$
is removed from $\mathit{Rel}$, then $a\in\mathrm{out}(C)$ (as $\delta(C)\cap
B\neq\emptyset$) and therefore also $a\in\mathrm{out}(D)$ (as $\mathit{Rel}$
was initialized to $\mathit{Rel}_{\mathit{I}\cap\mathit{Out}}$ on line 1 and
$(C,D)\in\mathit{Rel}$). Thus,
$D\cap(S\setminus\mathrm{\delta}^{-1}_{a}(S))=\emptyset$, which means that
$(C,D)$ is removed from $\mathit{Rel}$ even with the restricted
$\mathit{Remove}$. Therefore, it does not matter whether
$S\setminus\mathrm{\delta}^{-1}_{a}(S)$ is a subset of or it has an empty
intersection with $\mathit{Remove}$. ∎

As justified above, we can optimize LRT as follows. Sets
$\mathit{Remove}_{a}(B)$ are computed only if $a\in\mathrm{in}(B)$ and in that
case we only add states $q\in\mathrm{\delta}^{-1}_{a}(S)$ to
$\mathit{Remove}_{a}(B)$. As a result, we can reduce the number of required
counters by maintaining $\mathit{Count}_{a}(v,B)$ if and only if
$a\in\mathrm{in}(B)$ and $a\in\mathrm{out}(v)$.

##  5 Implementation and Complexity of OLRT

We first briefly describe the essential data structures (there are some
additional data structures required by our optimizations) and then we sketch
the complexity analysis. For the full details, see the technical report [9].

### Data Structures.

The input LTS is represented as a list of records about its states. The record
about each state $v\in S$ contains a list of nonempty
$\mathrm{\delta}^{-1}_{a}(v)$ sets111We use a list rather than an array having
an entry for each $a\in\Sigma$ in order to avoid a need to iterate over
alphabet symbols for which there is no transition., each of them encoded as a
list of its members. The partition $P$ is encoded as a doubly-linked list
(DLL) of blocks. Each block is represented as a DLL of (pointers to) states of
the block. Each block $B$ contains for each $a\in\Sigma$ a list of (pointers
on) states from $\mathit{Remove}_{a}(B)$. Each time when any set
$\mathit{Remove}_{a}(B)$ becomes nonempty, block $B$ is moved to the beginning
of the list of blocks. Choosing the block $B$ on line 5 then means just
scanning the head of the list of blocks. The relation $\mathit{Rel}$ is
encoded as a resizable boolean matrix.

Each block $B\in P$ and each state $v\in S$ contains an $\Sigma$-indexed array
containing a record $B.a$ and $v.a$, respectively. The record $B.a$ stores the
information whether $a\in\mathrm{in}(B)$ (we need the test on
$a\in\mathrm{in}(B)$ to take a constant time), If $a\in\mathrm{in}(B)$, then
$B.a$ also contains a reference to the set $\mathit{Remove}_{a}(B)$,
represented as a list of states (with a constant time addition), and a
reference to an array of counters $B.a.\mathit{Count}$ containing the counter
$\mathit{Count}_{a}(v,B)$ for each $v\in\mathrm{\delta}^{-1}_{a}(S)$. Note
that for two different symbols $a,b\in\Sigma$ and some $v\in S$, the counter
$\mathit{Count}_{a}(v,B)$ has different index in the array
$B.a.\mathit{Count}$ than the counter $\mathit{Count}_{b}(v,B)$ in
$B.b.\mathit{Count}$ (as the sets $\mathrm{\delta}^{-1}_{a}(S)$ and
$\mathrm{\delta}^{-1}_{b}(S)$ are different). Therefore, for each $v\in S$ and
$a\in\Sigma$, $v.a$ contains an index $v_{a}$ under which for each $B\in P$,
the counter $\mathit{Count}_{a}(v,B)$ can be found in the array
$B.a.\mathit{Count}$. Using the $\Sigma$-indexed arrays attached to symbols
and blocks, every counter can be found/updated in a constant time. For every
$v\in S,a\in\Sigma$, $v.a$ also stores a pointer to the list containing
$\mathrm{\delta}^{-1}_{a}(v)$ or $null$ if $\mathrm{\delta}^{-1}_{a}(v)$ is
empty. This allows the constant time testing whether $a\in\mathrm{in}(v)$ and
the constant time searching for the $\mathrm{\delta}^{-1}_{a}(v)$ list.

### Complexity analysis (Sketch).

We first point out how our optimizations influence complexity of the most
costly part of the code which is the main while loop. The analysis of lines
14–16 of LRT is based on the observation that for any two
$B^{\prime},D^{\prime}\in P_{\mathit{Sim}_{\mathit{I}}}$ and any $a\in\Sigma$,
it can happen at most once that $a$ and some $B$ with $B^{\prime}\subseteq B$
are chosen on line 14 and at the same time
$D^{\prime}\subseteq\mathit{Remove}_{a}(B)$. In one single iteration of the
while-loop, blocks $C$ are listed by traversing all
$\mathrm{\delta}^{-1}(v),v\in B$ (the $D$s can be enumerated during the
$\mathit{Split}$ operation). Within the whole computation, for any
$B^{\prime}\in P_{\mathit{Sim}_{\mathit{I}}}$, transitions leading to
$B^{\prime}$ are traversed on line 14 at most $P_{\mathit{Sim}_{\mathit{I}}}$
times, so the complexity of lines 14–16 of LRT is
$\mathcal{O}(\sum_{a\in\Sigma}\sum_{D\in P_{\mathit{Sim}}}\sum_{v\in
S}|\mathrm{\delta}^{-1}_{a}(v)|)=\mathcal{O}(|P_{\mathit{Sim}_{\mathit{I}}}||\delta|)$.
In the case of OLRT, the number and the content of remove sets is restricted
in such a way that for a nonempty set $\mathit{Remove}_{a}(B)$, it holds that
$a\in\mathrm{in}(B)$ and
$\mathit{Remove}_{a}(B)\subseteq\mathrm{\delta}^{-1}_{a}(S)$. Hence, for a
fixed $a$, $a$-transition leading to a block $B^{\prime}\in
P_{\mathit{Sim}_{\mathit{I}}}$ can be traversed only $|\\{D^{\prime}\in
P_{\mathit{Sim}_{\mathit{I}}}\mid a\in\mathrm{out}(D^{\prime})\\}|$ times and
the complexity of lines 14–16 decreases to $\mathcal{O}(\sum_{D\in
P_{\mathit{Sim}_{\mathit{I}}}}\sum_{a\in\mathrm{in}(D)}|\delta_{a}|)$.

The analysis of lines 17–20 of LRT is based on the fact that once $(C,D)$
appears on line 17, no $(C^{\prime},D^{\prime})$ with $C^{\prime}\subseteq
C,D^{\prime}\subseteq D$ can appear there again. For a fixed $(C,D)$, the time
spent on lines 17–20 is in $\mathcal{O}(\sum_{v\in
B}|\mathrm{\delta}^{-1}(v)|)$ and only those blocks $C,D$ can meet on line 17
such that $C\times D\subseteq\mathit{I}$. Thus, the overall time spent by LRT
on lines 17–20 is in $\mathcal{O}(\sum_{B\in
P_{\mathit{Sim}_{\mathit{I}}}}\sum_{v\in\mathit{I}(B)}|\mathrm{\delta}^{-1}(v)|)$.
In OLRT, blocks $C,D$ can meet on line 17 only if $C\times
D\subseteq\mathit{I}\cap\mathit{Out}$, and the complexity of lines 17–20 in
OLRT decreases to $\mathcal{O}(\sum_{B\in
P_{\mathit{Sim}_{\mathit{I}}}}\sum_{v\in(\mathit{I}\cap\mathit{Out})(B)}|\mathrm{\delta}^{-1}(v)|)$.

Additionally, OLRT refines $\langle
P_{\mathit{I}},\mathit{Rel}_{\mathit{I}}\rangle$ to $\langle
P_{\mathit{I}\cap\mathit{Out}},\mathit{Rel}_{\mathit{I}\cap\mathit{Out}}\rangle$
on line 1. This can be done by successive splitting according to the sets
$\mathrm{\delta}^{-1}_{a}(S),a\in\Sigma$ and after each split, breaking the
relation between blocks included in $\mathrm{\delta}^{-1}_{a}(S)$ and the ones
outside. This procedure takes time
$\mathcal{O}(|\Sigma||P_{\mathit{I}\cap\mathit{Out}}|^{2})$.

Apart from some other smaller differences, the implementation and the
complexity analysis of OLRT are analogous to the implementation and the
analysis of LRT [1]. The overall time complexity of OLRT is
$\mathcal{O}\bigl{(}|\Sigma||P_{\mathit{I}\cap\mathit{Out}}|^{2}+|\Sigma||S|+|P_{{\mathit{Sim}_{\mathit{I}}}}|^{2}+\sum_{B\in
P_{\mathit{Sim}_{\mathit{I}}}}(\sum_{a\in\mathrm{in}(B)}(|\mathrm{\delta}^{-1}_{a}(S)|+|\delta_{a}|)+\sum_{v\in(\mathit{I}\cap\mathit{Out})(B)}|\mathrm{\delta}^{-1}(v)|)\bigr{)}.$

The space complexity of OLRT is determined by the number of counters, the
contents of the $\mathit{Remove}$ sets, the size of the matrix encoding of
$\mathit{Rel}$, and the space needed for storing the $B.a$ and $v.a$ records
(for every block $B$, state $v$ and symbol $a$). Overall, it gives
$\mathcal{O}(|P_{\mathit{Sim}_{\mathit{I}}}|^{2}+|\Sigma||S|+\sum_{B\in
P_{\mathit{Sim}_{\mathit{I}}}}\sum_{a\in\mathrm{in}(B)}|\delta_{a}^{-1}(S)|)$.

Observe that the improvement of both time and space complexity of LRT is most
significant for systems with large alphabets and a high diversity of sets of
input and output symbols of states. Certain regular diversity of sets of input
and output symbols is an inherent property of LTSs that arise when we compute
simulations over tree automata. We address the impact of employing OLRT within
the procedures for computing tree automata simulation in the next section.

##  6 Tree Automata Simulations

In [1], authors propose methods for computing tree automata simulations via
translating problems of computing simulations over tree-automata to problems
of computing simulations over certain LTSs. In this section, we show how
replacing LRT by OLRT within these translation-based procedures decreases the
overall complexity of computing tree-automata simulations.

A (finite, bottom-up) _tree automaton_ (TA) is a quadruple
$A=(Q,\Sigma,\Delta,F)$ where $Q$ is a finite set of states, $F\subseteq Q$ is
a set of final states, $\Sigma$ is a ranked alphabet with a ranking function
$r:\Sigma\rightarrow\mathbb{N}$, and $\Delta\subseteq Q^{*}\times\Sigma\times
Q$ is a set of transition rules such that if $(q_{1}\ldots
q_{n},f,q)\in\Delta$, then $r(f)=n$. Finally, we denote by $r_{\mathsf{m}}$
the smallest $n\in\mathbb{N}$ such that $n\geq m$ for each $m\in\mathbb{N}$
such that there is some $(q_{1}\ldots q_{m},f,q)\in\Delta$. We omit the
definition of the semantics of TA as we will not need it, and we only refer to
[5, 1].

For the rest of this section, we fix a TA $A=(Q,\Sigma,\Delta,F)$. A _downward
simulation_ $D$ is a binary relation on $Q$ such that if $(q,r)\in D$, then
for all $(q_{1}\ldots q_{n},f,q)\in\Delta$, there exists $(r_{1}\ldots
r_{n},f,r)\in\Delta$ such that $(q_{i},r_{i})\in D$ for each $i:1\leq i\leq
n$. Given a downward simulation $D$ which is a preorder called an _inducing
simulation_ , an _upward simulation $U$ induced by $D$_ is a binary relation
on $Q$ such that if $(q,r)\in U$, then (i) for all $(q_{1}\ldots
q_{n},f,q^{\prime})\in\Delta$ with $q_{i}=q,1\leq i\leq n$, there exists
$(r_{1}\ldots r_{n},f,r^{\prime})\in\Delta$ with $r_{i}=r$,
$(q^{\prime},r^{\prime})\in U$, and $(q_{j},r_{j})\in D$ for each $j:1\leq
j\neq i\leq n$; (ii) $q\in F\implies r\in F$. From now on, let $\mathit{D}$
denote the maximal downward simulation on $A$ and $\mathit{U}$ the maximal
upward simulation on $A$ induced by $\mathit{D}$.

To define the translations from downward and upward simulation problems, we
need the following notions. Given a transition $t=(q_{1}\ldots
q_{n},f,q)\in\Delta$, $q_{1}\ldots q_{n}$ is its _left-hand side_ and
$t(i)\in(Q\cup\\{\square\\})^{*}\times\Sigma\times Q$ is an _environment_ —the
tuple which arises from $t$ by replacing state $q_{i}$, $1\leq i\leq n$, at
the $i^{th}$ position of the left-hand side of $t$ by the so called hole
$\square\not\in Q$. We use ${\it Lhs}$ of to denote the set of all left-hand
sides of $A$ and $Env$ to denote the set of all environments of $A$.

We translate the downward simulation problem on $A$ to the simulation problem
on the LTS
$A^{\bullet}=(Q^{\bullet},\Sigma^{\bullet},\\{\delta_{a}^{\bullet}\mid
a\in\Sigma^{\bullet}\\})$ where $Q^{\bullet}=\\{q^{\bullet}\mid q\in
Q\\}\cup\\{l^{\bullet}\mid l\in{\it Lhs}\\}$,
$\Sigma^{\bullet}=\Sigma\cup\\{1,\ldots,r_{\mathsf{m}}\\}\\}$, and for each
$(q_{1}\ldots q_{n},f,q)\in\Delta$, $(q^{\bullet},q_{1}\ldots
q_{n}^{\bullet})\in\delta^{\bullet}_{f}$ and $(q_{1}\ldots
q_{n}^{\bullet},q_{i}^{\bullet})\in\delta^{\bullet}_{i}$ for each $i:1\leq
i\leq n$. The initial relation is simply
$\mathit{I}^{\bullet}=Q^{\bullet}\times Q^{\bullet}$. The upward simulation
problem is then translated into a simulation problem on LTS
$A^{\odot}=(Q^{\odot},\Sigma^{\odot},\\{\delta^{\odot}_{a}\mid
a\in\Sigma^{\odot}\\})$, where $Q^{\odot}=\\{q^{\odot}\mid q\in
Q\\}\cup\\{e^{\odot}\mid e\in Env\\}$, $\Sigma^{\odot}=\Sigma^{\bullet}$, and
for each $t=(q_{1}\ldots q_{n},f,q)\in\Delta$, for each $1\leq i\leq n$,
$(q_{i}^{\odot},t(i)^{\odot})\in\delta^{\odot}_{i}$ and
$(t(i)^{\odot},q^{\odot})\in\delta^{\odot}_{a}$. The initial relation
$I^{\odot}\subseteq Q^{\odot}\times Q^{\odot}$ contains all the pairs
$(q^{\odot},r^{\odot})$ such that $q,r\in Q$ and $r\in F\implies q\in F$, and
$((q_{1}\ldots q_{n},f,q)(i)^{\odot},(r_{1}\ldots r_{n},f,r)(i)^{\odot})$ such
that $(q_{j},r_{j})\in\mathit{D}$ for all $j:1\leq i\neq j\leq n$. Let
$\mathit{Sim}^{\bullet}$ be the maximal simulation on $A^{\bullet}$ included
in $\mathit{I}^{\bullet}$ and let $\mathit{Sim}^{\odot}$ be the maximal
simulation on $A^{\odot}$ included in $\mathit{I}^{\odot}$. The following
theorem shows correctness of the translations.

######  Theorem 6.1 ([1])

For all $q,r\in Q$, we have
$(q^{\bullet},r^{\bullet})\in\mathit{Sim}^{\bullet}$ if and only if
$(q,r)\in\mathit{D}$ and $(q^{\odot},r^{\odot})\in\mathit{Sim}^{\odot}$ if and
only if $(q,r)\in\mathit{U}$.

The states of the LTSs ($A^{\bullet}$ as well as $A^{\odot}$) can be
classified into several classes according to the sets of input/output symbols.
Particularly, $Q^{\bullet}$ can be classified into the classes
$\\{q^{\bullet}\mid q\in Q\\}$ and for each $n:1\leq n\leq r_{\mathsf{m}}$,
{$q_{1}\ldots q_{n}^{\bullet}\mid q_{1}\ldots q_{n}\in{\it Lhs}\\}$, and
$Q^{\odot}$ can be classified into $\\{q^{\odot}\mid q\in Q\\}$ and for each
$a\in\Sigma$ and $i:1\leq i\leq r(a)$, $\\{t(i)^{\odot}\mid t=(q_{1}\ldots
q_{n},a,q)\in\Delta\\}$. This turns to a significant advantage when computing
simulations on $A^{\bullet}$ or on $A^{\odot}$ using OLRT instead of LRT.
Moreover, we now propose another small optimization, which is a specialized
procedure for computing $\langle
P_{\mathit{I}\cap\mathit{Out}}\mathit{Rel}_{\mathit{I}\cap\mathit{Out}}\rangle$
for the both of $A^{\odot}$, $A^{\bullet}$. It is based on the simple
observation that we need only a constant time (not a time proportional to the
size of the alphabet) to determine whether two left-hand sides or two
environments are related by the particular $\mathit{Out}$ (more specifically,
$(e_{1}^{\odot},e_{2}^{\odot})\in\mathit{Out}$ if and only if the inner
symbols of $e_{1}$ and $e_{2}$ are the same, and $(q_{1}\ldots
q_{n}^{\bullet},r_{1}\ldots r_{m}^{\bullet})\in\mathit{Out}$ if and only if
$n\leq m$).

### Complexity of the Optimized Algorithm.

We only point out the main differences between application of LRT [1] and OLRT
on the LTSs that arise from the translations described above. For
implementation details and full complexity analysis of the OLRT versions, see
the technical report [9].

To be able to express the complexity of running OLRT on $A^{\bullet}$ and
$A^{\odot}$, we extend $\mathit{D}$ to the set ${\it Lhs}$ such that
$((q_{1}\dotsc q_{n}),(r_{1}\dotsc r_{n}))\in\mathit{D}$ if and only if
$(q_{i},r_{i})\in\mathit{D}$ for each $i:1\leq i\leq n$, and we extend
$\mathit{U}$ to the set $Env$ such that
$((q_{1}\dotsc q_{n},f,q)(i),(r_{1}\dotsc
r_{n},f,r)(i))\in\mathit{U}{\iff}m=n\wedge
i=j\wedge(q,r)\in\mathit{U}\wedge(\forall k\in\\{1,...,n\\}.\ k\neq
i\Longrightarrow\ (q_{k},r_{k})\in\mathit{D})$. For a preorder $\rho$ over a
set $X$, we use ${X}/{\rho}$ to denote the partition of $X$ according to the
equivalence $\rho\cap\rho^{-1}$.

The procedures for computing $\mathit{Sim}^{\bullet}$ and
$\mathit{Sim}^{\odot}$ consist of (i) translating $A$ to the particular LTS
($A^{\bullet}$ or $A^{\odot}$) and computing the partition-relation pair
inducing the initial preorder ($\mathit{I}^{\bullet}$ or
$\mathit{I}^{\odot}$), and (ii) running a simulation algorithm (LRT or OLRT)
on it. Here, we analyze the impact of replacing LRT by OLRT on the complexity
of step (ii), which is the step with dominating complexity (as shown in [1]
and also by our experiments; step (ii) is much more computationally demanding
than step (i)).

As shown in the technical report [9], OLRT takes on $A^{\bullet}$ and
$\mathit{I}^{\bullet}$ space $\mathcal{O}(Space_{D})$ where
$Space_{D}{=}(r_{\mathsf{m}}+|\Sigma|)|{\it Lhs}\cup Q|+|{\it Lhs}\cup
Q/\mathit{D}|^{2}+|\Sigma||{\it
Lhs}/\mathit{D}||Q|+r_{\mathsf{m}}|Q/\mathit{D}||{\it Lhs}|$ and time
$\mathcal{O}(Space_{D}+|\Sigma||Q/\mathit{D}|^{2}+|{\it
Lhs}/\mathit{D}||\Delta|)$. On $A^{\odot}$ and $\mathit{I}^{\odot}$, OLRT runs
in time $O(Space_{U})$ where
$Space_{U}=(r_{\mathsf{m}}+|\Sigma|)|Env|+|Env/U|^{2}+|Env/U||Q|+|Q/U||Env|$
and time $\mathcal{O}(Space_{U}+|\Sigma||Q/U|^{2}+|Env/U||\delta|)$.

We compare the above results with [1], where LRT is used. LRT on $A^{\bullet}$
and $\mathit{I}^{\bullet}$ takes $\mathcal{O}(Space_{D}^{old})$ space where
$Space_{D}^{old}=(|\Sigma|+r_{\mathsf{m}})|Q\cup{\it Lhs}||{Q\cup{\it
Lhs}}/{\mathit{D}}|$, and $\mathcal{O}(Space_{D}^{old}+|\Delta||{Q\cup{\it
Lhs}}/{\mathit{D}}|)$ time. In the case of $A^{\odot}$ and $I^{\odot}$, we
obtain space complexity $\mathcal{O}(Space_{U}^{old})$ where
$Space_{U}^{old}=|\Sigma||Env||{Env}/{\mathit{U}}|$ and time complexity
$\mathcal{O}(Space_{U}^{old}+r_{\mathsf{m}}|\Delta||{Env}/{\mathit{U}}|)$.

The biggest difference is in the space complexity (decreasing the factors
$Space_{D}^{old}$ and $Space_{U}^{old}$). However, the time complexity is
better too, and our experiments show a significant improvement in space as
well as in time.

##  7 Experiments

We implemented the original and the improved version of the algorithm in a
uniform way in OCaml and experimentally compared their performance.

| LTS | LRT | OLRT
---|---|---|---
source | $|S|$ | $|\Sigma|$ | $|\delta|$ | time | space | time | space
random | 256 | 16 | 416 | 0.12 | 9.6 | 0.02 | 1.9
random | 4096 | 16 | 3280 | 13.82 | 714.2 | 2.02 | 78.2
random | 16384 | 16 | 26208 | o.o.m. | 268.85 | 4514.9
random | 4096 | 32 | 6560 | 62.09 | 1844.2 | 4.36 | 121.4
random | 4096 | 64 | 13120 | 158.38 | 3763.2 | 6.59 | 211.2
pc | 1251 | 43 | 49076 | 7.52 | 418.1 | 2.63 | 119.0
rw | 4694 | 11 | 20452 | 81.28 | 3471.8 | 19.25 | 989.3
lr | 6160 | 35 | 90808 | 390.91 | 12640.8 | 45.69 | 1533.6

Table 1: LTS simulation results

The simulation algorithms were benchmarked using LTSs obtained from the runs
of the abstract regular model checking (ARMC) (see [3, 4]) on several classic
examples—producer-consumer (pc), readers-writers (rw), and list reversal
(lr)—and using a set of tree automata obtained from the run of the abstract
regular tree model checking (ARTMC) (see [2]) on several operations, such as
list reversal, red-black tree balancing, etc. We also used several randomly
generated LTSs and tree automata.

| TA | LTS | LRT | OLRT
---|---|---|---|---
source | $|Q|$ | $|\Sigma|$ | $r_{\mathsf{m}}$ | $|\Delta|$ | $|S|$ | $|\Sigma|$ | $|\delta|$ | time | space | time | space
random | 16 | 16 | 2 | 245 | 184 | 18 | 570 | 0.06 | 6.2 | 0.02 | 1.4
random | 32 | 16 | 2 | 935 | 655 | 18 | 2165 | 0.87 | 74.4 | 0.21 | 14.4
random | 64 | 16 | 2 | 3725 | 2502 | 18 | 8568 | 26.63 | 1417.9 | 3.50 | 195.4
random | 32 | 32 | 2 | 1164 | 719 | 34 | 2511 | 2.67 | 166.6 | 0.23 | 16.8
random | 32 | 64 | 2 | 2026 | 925 | 66 | 3780 | 12.17 | 623.5 | 0.56 | 25.4
ARTMC111One of the automata selected from the ARTMC set. | 47 | 132 | 2 | 837 | 241 | 134 | 1223 | 0.84 | 70.6 | 0.05 | 6.2
ARTMC | variable222A set containing 10 305 tree automata of variable size (up to 50 states and up to 1000 transitions per automaton). The results show the total amount of time required for the computation and the peak size of allocated memory. | 517.98 | 116.2 | 80.84 | 22.1

Table 2: Downward simulation results

We performed the experiments on AMD Opteron 8389 2.90 GHz PC with 128 GiB of
memory (however we set the memory limit to approximately 20 GiB for each
process). The system was running Linux and OCaml 3.10.2.

The performance of the algorithms is compared in Table 1 (general LTSs), Table
2 (LTSs generated while computing the downward simulation), and Table 3 (LTSs
generated while computing the upward simulation), which contain the running
times ([s]) and the amount of memory ([MiB]) required to finish the
computation.

| TA | LTS | LRT | OLRT
---|---|---|---|---
source | $|Q|$ | $|\Sigma|$ | $r_{\mathsf{m}}$ | $|\Delta|$ | $|S|$ | $|\Sigma|$ | $|\delta|$ | time | space | time | space
random | 16 | 16 | 2 | 245 | 472 | 17 | 952 | 1.03 | 96.5 | 0.09 | 4.8
random | 32 | 16 | 2 | 935 | 1791 | 17 | 3700 | 18.73 | 1253.8 | 1.37 | 54.7
random | 64 | 16 | 2 | 3725 | 7126 | 17 | 14824 | 405.89 | 14173.9 | 22.83 | 752.6
random | 32 | 32 | 2 | 1164 | 2204 | 33 | 4548 | 64.10 | 3786.7 | 2.36 | 193.4
random | 32 | 64 | 2 | 2026 | 3787 | 65 | 7874 | o.o.m. | 6.72 | 245.8
ARTMC11footnotemark: 1 | 47 | 132 | 2 | 837 | 1095 | 133 | 3344 | 66.46 | 4183.2 | 0.69 | 68.2
ARTMC | variable22footnotemark: 2 | 12669.94 | 4412.6 | 400.62 | 106.6

Table 3: Upward simulation results

As seen from the results of our experiments, our optimized implementation
performs substantially better than the original. On average, it improves the
running time and space requirements by about one order of magnitude. As
expected, we can see the biggest improvements especially in the cases, where
we tested the impact of the growing size of the alphabet.

##  8 Conclusion

We proposed an optimized algorithm for computing simulations over LTSs, which
improves the asymptotic complexity in both space and time of the best
algorithm (LRT) known to date (see [1]) and which also performs significantly
better in practice. We also show how employing OLRT instead of LRT reduces the
complexity of the procedures for computing tree-automata simulations from [1].
As our future work, we want to develop further optimizations, which would
allow to handle even bigger LTSs and tree automata. One of the possibilities
is to replace existing data structures by a symbolic representation, for
example, by using BDDs.

Acknowledgements.  This work was supported in part by the Czech Science
Foundation (projects P103/10/0306, 102/09/H042, 201/09/P531), the BUT FIT
grant FIT-10-1, the Czech COST project OC10009 associated with the ESF COST
action IC0901, the Czech Ministry of Education by the project MSM 0021630528,
and the ESF project Games for Design and Verification.

## References

  * [1] P.A. Abdulla, A. Bouajjani, L. Holík, L. Kaati, and T. Vojnar.  Computing Simulations over Tree Automata: Efficient Techniques for Reducing Tree Automata.  In Proc. of TACAS’08, LNCS 4963. Springer, 2008.
  * [2] A. Bouajjani, P. Habermehl, L. Holík, T. Touili, and T. Vojnar.  Antichain-Based Universality and Inclusion Testing over Nondeterministic Finite Tree Automata.  In Proc. of CIAA’08, LNCS 5148. Springer, 2008.
  * [3] A. Bouajjani, P. Habermehl, P. Moro, and T. Vojnar.  Verifying Programs with Dynamic 1-Selector-Linked Structures in Regular Model Checking.  In Proc. of TACAS’05, LNCS 3440. Springer, 2005.
  * [4] A. Bouajjani, P. Habermehl, and T. Vojnar.  Abstract Regular Model Checking.  In Proc. of CAV’04, LNCS 3114. Springer, 2004.
  * [5] H. Comon, M. Dauchet, R. Gilleron, C. Löding, F. Jacquemard, D. Lugiez, S. Tison, and M. Tommasi.  Tree Automata Techniques and Applications.
http://www.grappa.univ-lille3.fr/tata, 2007.  release October, 12th 2007.

  * [6] D. Dams, O. Grumberg, and R. Gerth.  Generation of Reduced Models for Checking Fragments of CTL.  In Proc. of CAV’93, 1993.
  * [7] O. Grumberg and D. E. Long.  Model Checking and Modular Verification.  ACM Transactions on Programming Languages and Systems, 16, 1994\.
  * [8] M. R. Henzinger, T. A. Henzinger, and P. W. Kopke.  Computing simulations on finite and infinite graphs.  In Proc. of FOCS’95. IEEE Computer Society, 1995.
  * [9] L. Holík and J. Šimáček.  Optimizing an LTS-Simulation Algorithm.  Technical Report FIT-TR-2009-03, Brno University of Technology, 2009.
http://www.fit.vutbr.cz/~holik/pub/FIT-TR-2009-03.pdf.

  * [10] Lukáš Holík and Jiří Šimáček.  Optimizing an lts-simulation algorithm.  Computing and Informatics, 2010(7):1337–1348, 2010.
  * [11] R. Piage and R. Tarjan.  Three Partition Refinement Algorithms.  SIAM Journal on Computing, 16, 1987.
  * [12] F. Ranzato and F. Tapparo.  A New Efficient Simulation Equivalence Algorithm.  In Proc. of LICS’07, 2007.
  * [13] Z. Sawa and P. Jančar.  Behavioural Equivalences on Finite-State Systems are PTIME-hard.  Computing and Informatics, 24, 2005.