A Multi-Objective Multi-Modal Optimization Approach for Mining Stable Spatio-Temporal Patterns

1 Introduction

Spatio-temporal data mining is concerned with �nding speci�c patterns in databases describing temporally situated spatial objects. Many approaches have been developed in signal processing and computer science to address such a goal, ranging from Fourier Transforms to Independent Component Analysis [Hyvarinen et al., 2001], mixtures of models [Chudova et al., 2003] or string kernel machines [Saunders et al., 2004], to name a few. These approaches aim at particular pattern properties (e.g. independence, generativity) and/or focus on particular data characteristics (e.g. periodicity). However, in some application domains, the above properties are not relevant to extract the patterns of interest. The approach presented in this paper is motivated by such an application domain, functional brain imaging. Non invasive techniques, speci�cally magnetoencephalography (MEG) [H¤am¤al¤ainen et al., 1993], provide measures of the human brain activity with an unprecedented temporal resolution (time step 1 millisecond). This resolution enables researchers to investigate the timing of basic neural processes at the level of cell assemblies [Pantazis et al., 2005]. A key task is to �nd out which brain areas are active and when, i.e. to detect spatio-temporal regions with highly correlated activities. It is emphasized that such regions, referred to as stable spatio-temporal patterns (STPs), are neither periodic nor necessarily independent. Currently, STPs are manually detected, which (besides being a tedious task) severely hampers the reproducibility of the results. This paper addresses the automatic detection of stable spatio-temporal patterns, i.e. maximal spatiotemporal regions with maximal correlation. This detection problem cannot be formalized as a multi-objective optimization (MOO) problem [Deb, 2001], because experts are interested in several active brain areas: an STP might be relevant though it corresponds to a smaller spatio-temporal region, with a lesser correlated activity than another STP. The proposed approach thus extends MOO along the lines of multi-modal optimization [Li et al., 2002], de�ning a multiobjective multi-modal optimization framework (MoMOO). MoMOO is tackled by an evolutionary algorithm termed 4DMiner, using a diversity criterion to relax the multi-objective search (as opposed to diversity enforcing heuristics in MOO, e.g. [Corne et al., 2000; Laumanns et al., 2002]; more on this in section 2.3). To the best of our best knowledge, both the extension of evolutionary algorithms to MoMOO, and the use of multi-objective optimization within spatiotemporal data mining are new, although MOO attracts increasing attention in the domain of machine learning and knowledge discovery (see, e.g., [Ghosh and Nath, 2004; Francisci et al., 2003]). Experimental validation on arti�- cial and real-world datasets demonstrates the good scalability and performances of 4D-Miner; sought STPs are found within minutes from medium sized datasets (450 Mo), on PCPentium 2.4 GHz. The paper is organized as follows. Section 2 formalizes the detection of stable spatio-temporal patterns as a multi-modal multi-objective optimization problem, and motivates the use of evolutionary algorithms [B¤ack, 1995; Goldberg, 2002] for their detection. Section 3 gives an overview of 4D-Miner. Section 4 describes the experiment goals and setting. Section 5 reports on the extensive validation of 4D-Miner on arti�cial and real-world data. Section 6 discusses the approach with respect to relevant work, and the paper ends with perspectives for further research.

2 Position of the problem ; notations, criteria 2.1 Notations and de�nitions

Let N be the number of measure points. To the i-th measure point is attached a spatial position 1 Mi = (xi, yi, zi) and a temporal activity Ci(t), t ∈ [1, T]. Let I = [t1, t2] ⊂ [1, T]

1 MEG measure points actually belong to a 2D shape (the surface of the skull). However, the approach equally handles 2D or 3D spatio-temporal data.

be a time interval, and let ¯ C I i denote the average activity of the i-th measure point during I. The I-alignment σI (i, j) of measure points i and j over I is de�ned as:

As the sought patterns do not need to be spheric, ellipsoidal distances are considered. Only axis-parallel ellipsoids will be considered throughout the paper. For each weight vector w = (a, b, c) (w ∈ IR+3 ), distance dw is de�ned on the measure points as:

2.2 Multi-objective formalization

A pattern X = {i, I, w, r} is characterized by a center point i, i ∈ [1, N], a time interval I, an ellipsoidal distance dw, and a radius r > 0. The spatial neighborhoodB(i, w, r) is de�ned as the set of measure points j such that dw(i, j) is less than r.

The spatial amplitude of X, noted a(X), is the number of measure points in B(i,w, r). The temporal amplitude of X, noted ℓ(X), is the number of time steps in interval I. The spatio-temporal alignment of X, noted σ(X) is de�ned as the average, over all measure points in B(i, w, r), of their Ialignment with respect to the center i :

Interesting spatio-temporal patterns, according to the expert's �rst speci�cations, showmaximal spatial and temporal amplitudes together with a maximal spatio-temporal alignment. Speci�cally, a solution pattern is such that i) increasing its spatial or temporal amplitude would cause the spatiotemporal alignment to decrease; ii) the only way of increasing its spatio-temporal alignment is by decreasing its spatial or temporal amplitude. It thus comes naturally to formalize the STP detection in terms of multi-objective optimization problem (see [Deb, 2001] and references therein; [de la Iglesia et al., 2003; Bhattacharyya, 2000] for related data mining approaches).

The three a, ℓ and σ criteria induce a partial order on the set of patterns, referred to as Pareto domination.

De�nition 1. (Pareto domination)

Let c1, ...cK be K real-valued criteria de�ned on domain Ω, and let X and Y belong to Ω.

X Pareto-dominates Y , noted X ≻ Y , iff ck(X) is greater or equal ck(Y ) for all k = 1..K, and the inequality is strict for at least one k0.

The set of non-dominated solutions after a set of criteria is classically referred to as Pareto front [Deb, 2001].

2.3 Multi-modal multi-objective formalization

However, the above criteria only partially account for the expert's expectations: a STP might have a lesser spatiotemporal alignment and amplitude than another one, and still be worthy, provided that it corresponds to another active brain area. Therefore, not all sought STPs belong to the Pareto front.

Multi-modal optimization, interested in �nding global and local optima of the �tness landscape, has been extensively studied in the EC literature [Li et al., 2002]. Based on the multi-modal framework, a new framework extending multiobjective optimization and referred to as multi-modal multiobjective optimization (MoMOO), is thus proposed. Formally, let us �rst de�ne a relaxed inclusion relationship, noted p-inclusion.

De�nition 2. (p-inclusion)

Let A and B be two subsets of set Ω, and let p be a positive real number (p ∈ [0, 1]). A is p-included in B, noted

A ⊂p B, iff |A ⋂ B| > p× |A|. De�ning adequately the support of a candidate solution X (see below) multi-modal Pareto domination can be de�ned as follows:

De�nition 3. (multi-modal Pareto domination)

With same notations as in Def. 1, X mo-Pareto dominates Y , noted X ≻mo Y , iff the support of Y is p-included in that of X, and X Pareto-dominates Y .

In the case of STPs, the support of X = (i, I,w, r) is naturally de�ned as Supp(X) = B(i, dw, r) × I, with Supp(X) ⊂ [1, N] × [1, T ].

In contrast with [Corne et al., 2000; Laumanns et al., 2002] who use diversity-based heuristics for a better sampling of the Pareto front, diversity is thus used in MoMOO to rede�ne and extend the set of solutions.

2.4 Discussion

Functional brain imaging sets two speci�c requirements on spatio-temporal data mining. First, the expert's requirements are subject to change. Much background knowledge is involved in the manual extraction of stable spatio-temporal patterns, e.g. about the expected localization of the active brain areas. A �exible approach, accommodating further criteria and allowing the user to customize the search (in particular, tuning the thresholds on the minimal spatial or temporal amplitudes, or spatio-temporal alignment) is thus highly desirable.

Second, the approach must be scalable with respect to the data size (number of measure points and temporal resolution). Although the real data size is currently limited, the computational cost must be controllable in order to ef�ciently adjust the user-supplied parameters; in other words, the mining algorithm must be an any-time algorithm [Zilberstein, 1996].

The approach proposed is therefore based on evolutionary algorithms (EAs); these are widely known as stochastic, population-based optimization algorithms [B¤ack, 1995; Goldberg, 2002] that are highly �exible. In particular, EAs address multi-modal optimization [Li et al., 2002] and they can be harnessed to sample the whole Pareto front associated to a set of optimization criteria, with a moderate overhead cost compared to the standard approach (i.e., optimizing a weighted sum of the criteria gives a single point of the Pareto front) [Deb, 2001]. Last, the computational resources needed by EAs can be controlled in a straightforward way through limiting the number of generations and/or the population size.

3 Overview of 4D-Miner

This section describes the 4D-Miner algorithm designed for the detection of stable spatio-temporal patterns.

3.1 Initialization

Following [Daida, 1999], special care is devoted to the initialization of this evolutionary algorithm. The extremities of the Pareto front, where STPs display a high correlation (respectively, a low correlation) on a small spatio-temporal region (resp. a wide region), do not ful�ll the expert's expectations. Accordingly, some user-supplied thresholds are set on the minimal spatio-temporal amplitude and alignment.

In order to focus the search on relevant STPs, the initial population is generated using the initialization operator, sampling patterns X = (i, w, I, r) as follows:

• Center i is uniformly drawn in [1,N];
• Weight vector w is set to (1, 1, 1) (initial neighborhoods are based on Euclidean distance);
• Interval I = [t1, t2] is such that t1 is drawnwith uniform distribution in [1, T ]; the length t2-t1 of Ij is drawn according to a Gaussian distribution N(minℓ,minℓ/10), where minℓ is a time length threshold 2 .
• Radius r is deterministically computed from a usersupplied threshold minσ, corresponding to the minimal I-alignment desired.

All patterns X whose spatial amplitude a(X) is less than a user-supplied threshold mina are non-admissible; they will not be selected to generate new offspring. The user-supplied thresholds thus govern the proportion of usable individuals in the initial population.

2 In case t2 is greater than T, another interval I is sampled.

The computational complexity is in O(P × N ×minℓ), where P is the population size, N is the number of measure points andminℓ is the average length of the intervals.

3.2 Variation operators

From parent X = (i,w, I, r), mutation generates an offspring by one among the following operators:

• replacing center i with another measure point in B(i, w, r);
• mutating w and r using self-adaptive Gaussian mutation [B¤ack, 1995];
• incrementing or decrementing the bounds of interval I;
• generating a brand new individual (using the initialization operator).

The crossover operator, starting from parent X, �rst selects the mate pattern Y = (i′ , w′ , I ′ , r′ ) by tournament selection, Y minimizing the sum of the Euclidean distance between i and i′ , and the distance between the center of I and I ′ among K patterns uniformly selected in the population, where K is set to P/10 in all the experiments. The offspring is generated by:· replacing i with the center i′ of the mate pattern Y ;

• replacing w (resp. r) using an arithmetic crossover of w and w′ (respectively r and r′ );
• replacing I by the smallest interval containing I and I ′ .

An offspring is said admissible iff it satis�es the user-supplied thresholds mentioned in the initialization step.

3.3 Selection scheme

A Pareto archive is constructed with size L (set to 10 ×P in the experiments).

A steady-state scheme is used; at each step, an admissible parent X is selected among K uniformly drawn individuals in the population, by retaining the one which is dominated by the fewest individuals in the archive (Pareto tournament [Deb, 2001]). A single offspring Y is generated fromX by applying a variation operator among the ones mentioned above.

Offspring Y is evaluated by computing criteria a(Y ), ℓ(Y ), σ(Y ). Y is rejected if it is mo-Pareto dominated in the population; otherwise, it replaces a non-admissible individual in the population if any; if none it replaces an individual selected after anti-Pareto tournament (the individual out of K randomly selected ones in the population, that is dominated by the most individuals in the archive).

The archive is updated every P generations, replacing the mo-Pareto dominated individuals in the archive with individuals selected from the population after Pareto tournament.

4 Experimental setting and goal

This section presents the goal of the experiments, describes the arti�cial and real-world datasets used, and �nally gives the parameters of the algorithm and the performance measures used to evaluate the algorithm.

4.1 Goals of experiments and datasets

The initial goal is to provide the expert with a usable STP detection algorithm. The real datasets have been collected from people observing a moving ball. Each dataset involves 151 measure points and the number of time steps is 875. As can be noted from Fig. 1, which shows a representative dataset, the amplitude of the activities widely vary along the time dimension. The other goal is to assess the scalability and the performances of 4D-Miner, which is done using arti�cial datasets.

The arti�cial datasets are generated as follows. N measure points are located in uniformly selected locations in the 3D domain [0, 1]3. Activities are initialized from random cumulative uniform noise, with Ci(t + 1) = Ci(t) + ϵ, and ϵ is drawn according to U(0, τ). Every target STP S = (i,w, I, r, δ) is de�ned from a center i, a weight vector w, a time interval I, a radius r, and a fading factor δ, used to bias the activities as detailed below. The activity CS of the STP S is the average activity in the spatio-temporal region B(i, w, r) × I. Thereafter, activities are biased according to the target STPs: for each measure point j, for each STP S such that j is in�uenced by S (dw(i, j) < r), the activity Cj(t) is smoothed as

and d(t, I) is the distance of t to interval I (0 if t belongs to I, otherwise the minimum distance of t to the bounds of I). The scalability of 4D-Miner is assessed from its empirical

complexity depending on the number N of measure points and the number T of time steps.

The performances of 4D-Miner are measured using the standard criteria in information retrieval, recall and precision. The recall is the fraction of target STPs that are identi�ed by 4D-Miner, i.e. p-included in an individual in the archive; the precision is the fraction of relevant individuals in the archive (p-included in a target STP). For each experimental setting (see below), the recall and precision are averaged over 21 independent runs. The number of target STPs is set to 10 in all experiments. Each STP in�uences a number of measure points varying in [10,20], during intervals uniformly selected in [1, T ] with length varying in [15, 25], and δ uniformly selected in [0, .05]. Each problem instance thus includes STPs with various levels of dif�culty; the detection is hindered as the spatio-temporal support (r and I) is comparatively low and the δ parameter increases (the regularity is not visible outside interval I, as in Fig 2).

4.2 Experimental setting

The experiments reported in the next section considers a population size P = 200, evolved along 8000 generations (8,000 �tness evaluations per run). A few preliminary runs were used to adjust the operator rates; the mutation and crossover rates are respectively set to .7 and .3.

The number N of measure points ranges in {500, 1000, 2000, 4000}. The number T of time steps ranges in {1000, 2000, 4000, 8000}.

The thresholds used in the initialization are: mina = 5 (minimum number of curves supporting a pattern); minℓ = 5 (minimum temporal amplitude of a pattern); minσ = .1 (minimum spatio-temporal alignment of two curves in a pattern).

For computational ef�ciency, the p-inclusion is computed as: X is p-included in Y if the center i of X belongs to the spatial support of Y , and there is an overlap between their time intervals.

The maximal size of the datasets (T = 8000, N = 4000) is 456 Mo. Computational runtimes are given on PC-Pentium IV, 2.4 GHz. 4D-Miner is written in C++.

5 Experimental validation

This section reports on the experiments done using 4D-Miner.

5.1 Experiments on real datasets

Typical STPs found in the real datasets are shown in Figs. 3 and 4, showing all curves belonging to the STP plus the timewindow of the pattern. The runtime is less than 1 minute (PC Pentium 1.4 GHz). As discussed in section 2.3, many relevant STPs are Pareto-dominated: typically the STP shown in Fig 3 is dominated by the one in Fig 4.

These patterns are considered satisfactory by the expert. All experiments con�rm the importance of the user-de�ned thresholds, controlling the quality of the initial population. Due to the variability of the data, these threshold parameters are adjusted for each new experiment.

The coarse tuning of the parameters can be achieved based on the desired proportion of admissible individuals in the initial population. However, the �ne-tuning of the parameters could not be automatized up to now, and it still requires running 4D-Miner a few times. For this reason, the control of the computational cost through the population size and number of generations is one of the key features of the algorithm.

5.2 Multi-objective multi-modal optimization

The good scalability of 4D-Miner is illustrated in Fig. 5. The empirical complexity of the approach is insensitive to the number of time steps T and linear in the number N of measure points. This computational robustness con�rms the analysis (section 3.1), the evaluation of each pattern has linear complexity in N. On-going work is concerned with exploiting additional data structures inspired from [Wu et al., 2004] in order to decrease the complexity in N. Table 1 reports the

recall achieved by 4D-Miner over the range of experiments, averaged over 21 independent runs, together with the standard deviation. On-line performances are illustrated on Fig

4000 89 ± 10 81 ± 13 56 ± 14 32 ± 16

Table 1: Recall achieved by 4D-Miner vs N and T (average percentage and standard deviation over 21 runs)

6. These results con�rm the robustness of the proposed approach: a graceful degradation of the recall is observed as T and N increase. It must be noted that STPs occupy a negligible fraction of the spatio-temporal domain (circa 10-4 for T = 8000, N = 4000). The average precision is low, ranging from 12 to 20% (results omitted due to space limitations). However, post-pruning can be used to sort the wheat from the chaff in the �nal archive, and increase the precision up to

50% without decreasing the recall; the post-pruning straightforwardly removes the STPs with small spatial or temporal amplitudes.

Quite different effects are obtained when the archive is pruned along the search (e.g. by increasing the thresholds on the minimal spatial and temporal amplitudes), which decreases the overall performances by an order of magnitude; interestingly, similar effects are observed in constrained evolutionary optimization when the fraction of admissible solutions is very

low. A �nal remark is that the performances heavily depend upon the user-supplied thresholds (section 3.1), controlling the diversity and effective size of the population. Indeed, a parametric model of the dataset would enable automatic adjustment of these parameters. It might also be viewed as a advantage of 4D-Miner that it does not require any prior model of the data, and that inexpensive preliminary runs can be used to adjust the needed parameters.

6 Relevant work

This brief review does not claim exhaustiveness; the interested reader is referred to [Shekhar et al., 2003; Roddick and Spiliopoulou, 2002] for comprehensive surveys. Spatiotemporal data mining applications (e.g., remote sensing, environmental studies, medical imaging) are known to be computationally heavy; as most standard statistical methods do not scale up properly, new techniques have been developed, including randomized variants of statistical algorithms. Many developments are targeted at ef�cient access primitives and/or complex data structures (see, e.g., [Wu et al., 2004]); another line of research is based on visual and interactive data mining (see, e.g., [Keim et al., 2004]), exploiting the unrivaled capacities of human eyes for spotting regularities in 2Ddata. Spatio-temporal data-mining mostly focuses on clustering, outlier detection, denoising, and trend analysis. For instance, [Chudova et al., 2003] used EM algorithms for non-parametric characterization of functional data (e.g. cyclone trajectories), with special care regarding the invariance of the models with respect to temporal translations. The main limitation of such non-parametric models, including Markov Random Fields, is their computational complexity, sidestepped by using randomized search for model estimates.

7 Discussion and Perspectives

This paper has proposed a stochastic approach for mining stable spatio-temporal patterns. Indeed, a very simple alternative would be to discretize the spatio-temporal domain and compute the correlation of the signals in each cell of the discretization grid. However, it is believed that the proposed approach presents several advantages compared to the brute force, discretization-based, alternative. First of all, 4D-Miner is a fast and frugal algorithm; its good performances and scalability have been successfully demonstrated on medium sized arti�cial datasets. Second, data mining applications speci�- cally involve two key steps, exempli�ed in this paper: i) understanding the expert's goals and requirements; ii) tuning the parameters involved in the speci�cations. With regard to both steps, the ability of working under bounded resources is a very signi�cant advantage; any-time algorithms allow the user to check whether the process can deliver useful results, at a low cost. Further research is concerned with extending 4D-Miner in a supervised learning perspective (�nding the STPs that are complete - active in several persons undergoing the same experiment - and correct - not active in the control experiment). The challenge is to directly handle the additional constraints of completeness and correction in the multi-objective multi-modal optimization framework presented here.

Acknowledgments

The authors are partially supported by the PASCAL Network of Excellence (IST-2002-506778); this study was supported by a grant from the French Ministry of Research & Technology �ACI, New Applications of Mathematics�.