Parallelizing Convergent Cross Mapping Using Apache Spark

TL;DR

This paper presents a method to accelerate convergent cross mapping (CCM) using Apache Spark, enabling faster causal analysis in complex systems with large time series data.

Contribution

It introduces parallelized solutions for CCM on Spark, achieving significant speedups and scalability for causal inference in complex systems.

Findings

Achieved over tenfold performance improvement with parallelized CCM.

Demonstrated high scalability of the Spark-based CCM implementation.

Enabled faster causal analysis in large, complex datasets.

Abstract

Identifying the causal relationships between subjects or variables remains an important problem across various scientific fields. This is particularly important but challenging in complex systems, such as those involving human behavior, sociotechnical contexts, and natural ecosystems. By exploiting state space reconstruction via lagged embedding of time series, convergent cross mapping (CCM) serves as an important method for addressing this problem. While powerful, CCM is computationally costly; moreover, CCM results are highly sensitive to several parameter values. While best practice entails exploring a range of parameter settings when assessing casual relationships, the resulting computational burden can raise barriers to practical use, especially for long time series exhibiting weak causal linkages. We demonstrate here several means of accelerating CCM by harnessing the distributed Apache Spark platform. We characterize and report on results of several experiments with parallelized solutions that demonstrate high scalability and a capacity for over an order of magnitude performance improvement for the baseline configuration. Such economies in computation time can speed learning and robust identification of causal drivers in complex systems.

Tables

Table 1. Table 1: Implementation Levels

	Implementation Level
Case A1	Single-threaded CCM (no RDD & Pipeline)
Case A2	Synchronous CCM Transform Pipelines
Case A3	Asynchronous CCM Transform Pipelines
Case A4	Synchronous Distance Indexing Table & CCM Transform Pipelines
Case A5	Asynchronous Distance Indexing Table & CCM Transform Pipelines

Table 2. Table 2: Elasticity Analysis

Parameter varied	parameter	Case B1	Case B2	Case B3
$L$	$L$	500	1000	2000
	others	the same as baseline scenario
$E$	$E$	1	2	4
	others	the same as baseline scenario
$\tau$	$\tau$	1	2	4
	others	the same as baseline scenario

Full text

1 Introduction

Identification of causal relations between variables in many domains has traditionally relied upon controlled experimentation, or investigation of underlying mechanisms. The first of these requires a heavy investment of time and financial resources, and can pose ethical concerns. The limits of such controlled studies – such as Randomized Controlled Trials (RCTs) – are particularly notable in the context of complex systems, which exhibit reciprocal feedbacks, delays and non-linearities [6]. While ubiquitous in science, progress in the latter approach is commonly measured in decades. In recent years, a growing number of researchers have applied Pearl’s causal inference framework [11] to identify causal linkages, but such methods can be challenging to apply in the presence of observability constraints, and when the variables of interest are coupled but distant within the system. Most critically, in the context of complex systems, such inference techniques encounter challenges in the context of reciprocal causality.

Convergent cross mapping [13] is an algorithm based on Takens’ embedding theorem [14] that can detect and help quantify the relative strength of unidirectional and bidirectional causal relationships between variables X and Y in coupled complex systems. Within the CCM algorithm, in order to assess if variable $Y$ is causally governed by variable $X$, we attempt to predict the value of $X$ on the basis of the state space reconstructed from $Y$ (see below); for statistical reliability, this must be done over a large number $r$ of realizations. To assess causality, we examine whether these results ”converge” as we consider a growing numbers $L$ of datapoints within $Y$ within our reconstruction (also see below).

Overall, the results of CCM are sensitive to the parameters below:

• $r$ The number of random subsamples, commonly 250 or larger.

• $\tau$ The embedding delay used in the shadow manifold reconstruction.

• $E$ The embedding dimension of the dynamic system. For simplex projection, $E$ will range from 1 to 10.

• $L$ The size of the library extracted from time series.

Running CCM with a wide range of different parameter settings imposes a high computational overhead. But, as for many data science tasks, we believed that the performance could be elevated via parallel and distributed processing by implementing the CCM algorithm atop Apache Spark [19] (henceforth, “Spark”) and distributing computations across a Yarn cluster [15].

In this paper, following additional background on CCM and the related literature, we describe a CCM parallel implementation which utilizes the MapReduce framework [2] provided by Spark. The paper then presents a performance evaluation and comparison of the framework. We can conclude from the experiments that, with the parallel techniques and cloud computing support, researchers can use CCM to confidently infer causal connections between larger time series in far less time than previously required.

2 Background

2.1 Convergent Cross Mapping Basics

In 2012, Sugihara et al. [13] built on ideas from Takens’ Theorem [14] to propose convergent cross mapping (CCM) to test causal linkages between nonlinear time series observations. This approach has enjoyed varied applications. For example, Luo et al. [7] successfully revealed underlying causal structure in social media and Verma et al. [16] studied cardiovascular and postural systems by taking advantages of this algorithm.

We provide here a brief intuition for why and how CCM works. Consider two variables $X$ and $Y$, each associated with eponymous time-series and – further – where $Y$ depends on $X$. For example, consider a case where for each timepoint $X$ measures the count of hares, and $Y$ that of lynx. In this situation, if $Y$ (lynx) causally depends on $X$ (hares), observing the values of $Y$ over time (e.g., a steep drop or a plateauing in lynx numbers) tells us about the state of governing factors, including $X$ (here, the fact that the number of hares is too small to effectively feed the lynx population, or that they are roughly in balance with lynx, respectively). A implication of this – captured by Takens’ Theorem – is that information on the state of $X$ is encoded in the state space reconstructed from $Y$, meaning that points that are located nearby within $Y$’s reconstructed state space will be associated with similar values for $X$, and can thus be used to make accurate (skillful) prediction of the value of $X$. In most cases, such prediction of one variable (e.g., $X$) within the state space of another ($Y$) can be achieved by nearest neighbor forecasting using simplex projection [3]. Pearson correlation between observed and predicted values can be applied to measure prediction skill.

2.2 Past work in CCM Performance Improvement

Despite the fact that CCM is increasingly widely applied, there remain pronounced computational challenges in applying the tool for moderate and large time series. In order to secure confidence in inferences regarding causality, use of appropriate parameter values and a relative longer input are required for the original CCM [10]. As such, since its first appearance in 2012, a number of modifications and improvements have been proposed to handle this drawback. In 2014, Ma et al. [8] developed cross-map smoothness (CMS) based on CCM which has the advantage of requiring a shorter time series. Compared to original CCM, CMS can be used for time series in the order of $n=10$, whereas CCM requires time-series at least in the order of $n=10^{3}$ to yield reliable results. Additionally, works such as [1], [4], [5] investigated and introduced mathematical methods to properly estimate parameters required by CCM (embedding dimension $E$, time delay $\tau$ and time subsample $L$). Such work expanded CCM-related research and also provided methods for quickly inferring causality in certain circumstances.

The previous improvements on CCM typically trade off potential accuracy for relatively fast execution, and the assumptions in some methods cannot be safely maintained with noisy time series observations. However, the original CCM can be improved by introduction of parallel computing techniques. In recent years, numerous studies such as [9], [12] have been conducted using distributed computing frameworks such as MPI or Spark. Such parallel techniques can dramatically improve the algorithmic performance by effectively exploiting the cluster-based computational capacity. It is worthwhile to implement a parallel version of CCM to allow researchers to rapidly and robustly evaluate the existence and strength of causal connections between measured time series.

3 Methodology

To achieve a Spark parallel version of CCM, we introduce two core concepts: the Spark Resilient Distributed Dataset (RDD) [18] and Pipeline. The former is the immutable data structure that can be operated in a distributed manner, which brings significant benefits for concurrently draw $r$ subsamples of time series to assess Cross-Mapping convergence. As for the pipeline, it is specified as a sequence of stages, and each stage transforms the original RDD to another RDD accordingly. In summary, the definition of pipeline supports an elegant design for a parallel CCM algorithm manipulating RDDs in Spark.

3.1 CCM Transform Pipeline

Consider applying CCM to test if the variable associated with time series $Y$ is being driven by the variable associated with time series $X$. In the corresponding transform pipeline, the parallel version of CCM is implemented as several stages to transform the RDD of $r$ random subsamples of the time series to the RDD of prediction skills for a given ($\tau$, $E$, $L$) tuple. To start the transformation, an input RDD is created which includes a pair of subsamples of lengths $L$ of each of the time series, and values for each of two parameters ($\tau$, $E$). The output of the CCM transform pipeline is an RDD of sequences of prediction skills. In the whole procedure, Spark operates the whole transformation in parallel without extra coding as shown in Fig. 1.

3.2 Distance Indexing Table Pipeline

The CCM transform pipeline above achieves the aim of running CCM concurrently on multiple subsamples $r$. However, there is still a considerable potential for further optimization for this pipeline. As mentioned before, the most time-consuming part in the original CCM is finding the $E+1$ nearest neighbors for every lagged-coordinate vector ($\tau$) in the shadow manifold. For every point in the input RDD, the CCM transform pipeline computes the distances to all lagged-coordinate vectors of subsamples, sorts them and finally takes the top $E+1$ as the nearest neighbors. This process is inefficient because of its repeated sorting and calculation for all the subsamples. It is particularly notable that as the length of subsamples $L$ used for computation increases, the running time will grow superlinearly.

The best approach is to break down the nearest neighbors searching of CCM transformations into two parts: distance indexing table construction and nearest neighbors searching based on the constructed table. The first part can be achieved by setting another pipeline as a preprocessing step before applying the CCM transform pipeline. After building the distance indexing table, Spark can broadcast this table to each worker node on the cluster at one time rather than ship a copy of it every time they need it, as shown in Fig. 2. The pipeline of constructing the distance indexing table will be executed concurrently on the entire input time series, and it also reduces a significant amount of repeated calculation in the CCM transform pipeline. From the experiment results, the total computation time decreases in a pronounced fashion. As the library size $L$ grows, the time spent on searching for the nearest neighbors increases correspondingly, and pre-building the distance indexing table secures increasing benefit.

3.3 Asynchronous Pipelines

After the pipeline is created to run CCM, a job is generated in the master node and then submitted to the cluster and partitioned into many tasks running in the executors of worker nodes. This setting is defined in the job submission and is, in general, constant. If we perform two pipelines one after other, they will always be executed sequentially. As such, we can adopt the asynchronous mechanisms to increase the parallelism and execute different pipelines concurrently. FutureAction is the Spark API to undertake asynchronous job submission. It provides a native way for the program to express concurrent pipelines without having to deal with the detailed complexity of explicitly setting up multiple threads. In this way, we can achieve running various combinations of the parameters ($L$, $\tau$, and $E$) in parallel by executing multiple concurrent pipelines.

4 Experiment Results

The baseline scenario of parameters, with input time series size of 4000, $r$ of 500, $L$ with values $[500,1000,2000]$, $E$ and $\tau$ both with $[1,2,4]$, is set for the comparison in the experiments. In the following experiments, the Spark parallel version of CCM will be run three times on the Google Cloud Platform to obtain the average computation time. The cluster setting is 1 master node and 5 worker nodes with 4 cores CPU and 15 GB Memory.

4.1 Overview of Improvements

This experiment compares the performance improvement of different implementations on the baseline scenario. These implementations in Table 1 are submitted on Yarn Mode and Local Mode, separately. Yarn Mode, or cluster mode will exploit all the worker nodes existing in the cluster while Local Mode only runs applications on the master node (Single Machine).

The results are shown in Fig 4. Several conclusions can be drawn from the experimental results for different levels of parallel implementation. Firstly, the single-threaded version of CCM imposes a heavy computational cost, and there is no difference between two modes as they do not utilize the worker nodes in the cluster. Next, asynchronous pipelines can only reduce computation time in Yarn mode. After the comparison of the CPU utilization rates, it indicates that the asynchronous pipelines cannot offer more parallelization when the CPU utilization already reaches full throttle. However, when run with Yarn, the worker nodes still have rooms to improve the utilization rates. Also, as seen from the results, the spark full parallel version (Case A5) is approximately 1.2% the running time of the single-threaded version. Ultimately, the most significant improvement of marginal computation performance lies in adding the distance indexing table pipeline based on the CCM Transform pipeline. It reduces the computation time cost by over 80% relative to the baseline. Such considerable improvement shows the parallel version of CCM benefits strongly from establishing the distance indexing table globally for nearest neighbors searching pipeline.

When comparing current existing public CCM implementation, rEDM R package, which created by the Hao Ye et al. [17] using lower level language C++, our Spark parallel implementation (Case A5) is approximately 15x faster than rEDM for baseline scenario on current cluster setup on Google Compute Platform. Obviously, the parallel version can perform more favorably with a more powerful cluster (vertical scaling) or adding more workers (horizontal scaling).

4.2 Elasticity Analysis

As a range of parameter settings been looped over for the best results to infer causality, testing the elasticity of running time concerning a given parameter value is necessary. Two versions of CCM (Parallel Version is the implementation Case A5 which has all degrees of parallelism, while Single-threaded Version is the implementation Case A1 which has not implemented any pipeline) will be tested with the parameter settings as shown in Table 2. Intuitively, these cases vary only one parameter from the baseline for comparison. When doubling parameter $L$, the average run time increases to 4.06x using the Spark single-threaded version, and 1.11x using the Spark parallel version. Similarly, doubling parameter $\tau$ and $E$ almost has no impact on running time in the parallel version. However, doubling $\tau$ indeed increases the running time to 1.13x in the single-threaded version as it increase the dimension of shadow manifolds $M_{x}$ and $M_{y}$.

In summary, the values of these parameters, especially for $L$, do influence execution time for both the single-threaded and parallel versions; however, with the current optimization of the parallel methods by introducing indexing table before nearest neighbor searching, the relative impact shrinks, which make testing relatively large parameters for the causality assessment a reality.

5 Conclusion

The Spark framework provides relatively convenient APIs to exploit parallelism in algorithms such as CCM. This work conducted experiments demonstrating the performance benefits of exploiting the parallelism in CCM algorithm using Spark. The scalability of Spark offers considerable benefits in accelerating the execution with the support of clusters, allowing for a significant reduction in running time when adding more worker nodes into the cluster. Of critical importance for robust application of Spark, these performance gains make this algorithm a valuable modeling tool to assess causality with confidence in an abbreviated time. Such gains are particularly important in the context of high-velocity datasets involving human behavior and exposures, such as are commonly collected in human social and sociotechnical systems.

While it demonstrated potential for marked speedups, this work suffers from some pronounced limitations. Construction of the distance indexing table trades off higher space consumption for savings in computation time; for large shadow manifolds from a large value of $L$, the indexing table can require a large quantities of system memory. However, as previous study [8] shows, CCM can produce reliable results when input time-series length is around in the order of $n=10^{3}$, for which the required memory space is well under what most current hardware can offer.