BACKGROUND
the exponential increase and ubiquitous use of computational technology has given rise to an era of "big data" that pushes the limits of conventional data analysis  <cit> .
techniques for analyzing big datasets often proceed by identifying patterns of co-occurrence or correlation through principal component analysis   <cit> , multidimensional scaling   <cit> , etc.
however, many of these methods require significant data reduction or smoothing which makes them difficult to interpret  <cit> .
other methods such as multiple linear regression or pearson's correlation coefficient  are easy to interpret as they operate on data in their native data space, without any kind of large data transformation or dimensionality reduction, but are limited in the structure that they can detect.
though pcc is a classic and powerful technique for finding linear relationships between two variables, it is not designed for capturing lead-lag relationships seen in time series data.
local similarity analysis   <cit>  extends correlation calculations to include the time variable, enabling identification of local correlates.
furthermore, ruan et al.
have presented a graphical network framework in which to visualize the structure of significant lsa correlations.
unfortunately, the current implementation of lsa requires multiple runs on permuted data and a monte carlo statistical method known as a permutation test to evaluate a null distribution and obtain a p-value determining significance.
each iteration of this procedure has a computational complexity of o, where p is the number of permutations, m is the number of covariate time series, and n is their length.
due to the number of pair-wise calculations needed, extant lsa is computationally onerous when m is large, limiting its use to datasets where the number of observed variables at each time point is small .
though there has been some improvement to its performance  <cit> , assumptions of normality and implementation issues continue to stymie practical application of lsa on big datasets.
here we describe a novel asymptotic upper bound on the calculation of the lsa statistic's p-value, resulting in an exponentially converging calculation to bound and check for significance of computed lsa statistics without a computationally intensive permutation test.
this bound does not require a rank-vector normal transformation, promoting application to any distribution that has finite variance.
as a result, this implementation of lsa can navigate big datasets with millions of co-variate time series.
we demonstrate this using time-series datasets from public health  <cit> , microbial ecology  <cit> , and social media  <cit> .
the implemented algorithm, named fastlsa, is written in c and optimized for threading on multi-core computers.
interpreting the lsa statistic
lsa concerns itself with pairs of time series data.
the lsa statistic can be interpreted in a manner similar to pcc when no lag window exists between two time series.
however, lsa is also capable of capturing localized correlation that is staggered or lagged.
a large positive or negative lsa value indicates a correspondingly strong pcc correlation or a correlation at a time displacement within the lag window .
note that if both positive and negative correlations exist between two time series, lsa will only report the strongest of the two.
lsa is advantageous on large datasets containing many time series.
results can be visualized as a graphical network where nodes represent the individual time series and the edges represent their lsa correlation statistic.
when displayed using a force-directed layout in cytoscape  <cit> , closely related time series cluster together, visually isolating clusters of local similarity.
metadata related to experimental or environmental conditions can then be applied to the nodes, shedding insight into hierarchical network structure.
implementation
description of the lsa algorithm
in this section we reproduce the algorithm from  <cit>  to compute lsa statistics and their corresponding p-values between pairs of time series in a dataset.
we assume as input a set of time series vectors of equal length.
let us denote the number of time series by m and their length as n. let us denote the time series dataset as x where xij denotes the jth element of the ith time series, with i =  <dig>   <dig>  ..., m and j =  <dig>   <dig>  ..., n, and assume that the xij are real numbers.
we also assume that there are no missing values in the dataset x, and realize that practical use will require interpolation or filtering.
in figure  <dig> we present the algorithm for computing the lsa statistic for a pair of time series, x=xi1n and y=yi1n, where the length of the time series is assumed to be equally spaced in time.
we have modified the presentation of the lsa algorithm by  <cit>  to highlight our analysis and derivation of a bound on the tail distribution of the lsa statistic.
specifically we calculate the lsa statistic for a pair of time series, x and y. two-dimensional arrays pi, j and ni, j are used to store the positive and negative partial sums  of the pairwise product of time series values.
we also assume a suitable bound on the maximum time lag considered while computing the lsa statistic, denoted by d.
the algorithm first initializes the arrays pj, <dig>  nj,  <dig>  p <dig>  i, and n <dig> i for all i, j =  <dig>  ..., n, with a maximum absolute difference of d. next it considers the time series pairs for each possible lag, up to a maximum of d, and then computes the progressive sum of the pair-wise products of the time-series values from the low to high index of the arrays.
during the computation, the progression of the partial sum is reset to  <dig> if the sum is below  <dig>  after partial sums have been computed, the values of n ^ and p ^ are calculated by taking the maximum of the corresponding values of the arrays n and p. finally, the lsa statistic is estimated as signp ^-n ^maxp ^,n ^n.
calculating the upper bound
in this section we derive the asymptotic upper bound on the p-value for the cumulative probability distribution of the lsa statistic without the need of a normality assumption.
our derivation is based on distributional results of the maximum cumulative sum of independent random variables known in the literature from probability theory  <cit> .
we begin by stating our assumptions about the dataset, isolate target calculations from the lsa algorithm, and from our referenced mathematical results, derive and prove important lemmas.
these lemmas will serve as the building blocks as we logically construct a theorem which will form the basis of our lsa p-value upper bound.
we begin by making certain assumptions about the probability model used to derive the bounds.
first, each pi, j or ni, j is considered individually.
we assume that the time series values xi, yj for i, j =  <dig>  ..., n are independent of one another.
this assumption can be made when weak dependence exists because it is near the truth and effective, much like the naive bayes assumption.
this assumption is also enabling, as it allows us to invoke the distributions of partial sums of independent random variables and continue in a mathematically straightforward way.
further, we assume independence between each time time series as a null hypothesis, and as it is subject to rejection upon obtaining a statistically significant lsa value.
consider lines  <dig> and  <dig> of the lsa algorithm , pi+k+ <dig> j+k+ <dig> ← max{ <dig>  pi+k, j+k + xi+k * yj+k} and p ^←maxi,j:|i-j|≤dpi,j.
for any pair of i and j let us define the sequence random variables as zk = xi+kyj+k for k =  <dig>  ..., min{n - i, n - j} -  <dig>  and the sequence of random variables ζk = z <dig> + ... + zk for k =  <dig>  ..., min{n - i, n - j} -  <dig> supposing ζ <dig> =  <dig>  using the above ζk's, we define random variables ηk* as ηk*=maxζ <dig> ζ <dig> ⋯,ζk for the same values of k =  <dig>  ..., min{n - i, n - j} -  <dig>
we also define the set of random variables η <dig>  η <dig>  ..., ηk by the recurrence formula ηk+ <dig> = max{ <dig>  ηk + zk+1}.
note that the random variables pi+k, j+k and the ηk have the same distribution.
it is shown in  <cit>  that the random variables ηk* and ηk also have the same distribution.
as a result, now we can analyze the cumulative distribution of pi+k, j+k as a distribution for ηk*, and use the results by nevzorov and petrov  <cit>  on pi+k, j+k to derive tail probability bounds.
we also assume that the random variables zk have the first two moments, although such assumptions are not required for the results of  <cit> , we use them to derive simpler bounds.
we now consider a few useful lemmas that we will use to construct our p-value upper bound.
the first step is to simplify the tail event  into simpler terms.
the following lemma expresses the tail event for lsa {|lsa| >x} and any x∈ℝ in terms of the tail events of {pi, j >x} and {ni, j >x}, the positive and negative lsa calculations for the same x, the bound on our test statistic .
lemma  <dig> for any x∈ℝwe have {|lsa| >x} = { ∪ }.
proof.
the result is clear from the following:
|lsa|>x=maxp ^,n ^>xn=p ^≤xn∩n ^≤xnc=maxijpij≤xn∩maxijnij≤xnc=∩ijpij≤xn∩∩ijnij≤xnc=∪ijpij>xn∪∪ijnij>xn □
in the lsa algorithm, we have maximums pij = max{ <dig>  pi- <dig>  j- <dig> + xi-1yj-1} and nij = max{ <dig>  pi- <dig> j- <dig> - xi- <dig>  j-1}, which complicates their theoretical analysis.
fortunately, equivalence have been demonstrated in the literature  <cit> , and we restate these in the following lemma for clarity: the similarity of the distributions of ηk and ηk*, for k =  <dig>  ..., min{n - i, n - j} -  <dig>  this will help us derive the bounds for the events {pij >xn} and {nij >xn}, the simpler terms we derived in the previous lemma.
lemma  <dig> let zi be mutually independent random variables and let us denote by sk= ∑i=1kziwhere s <dig> =  <dig>  and qk+ <dig> = max{ <dig>  qk + zk} with q <dig> =  <dig>  then p = p  for x∈ℝ.
in order to get a simple formula for the bound on the cumulative tail probabilities for pi, j and ni, j we reproduce below the results on partial sums of random variables due to nevzorov and petrov  <cit> .
for our sequence of independent and identically distributed  random variables under consideration {xn} it follows that lindeberg's condition holds  <cit> .
a property showing the variance of a distribution stabilizes as more variables are added, pinning the tails of it down.
thinking about this in terms of time series, as a series gets larger, the upper bound of the distribution becomes more defined and calculable.
now to build theorems upon which we will derive a formulaic p-value bound.
theorem  <dig> if the random variables {xn} have zero expectation and finite variances and if lindeberg's condition holds: Λn →  <dig> as n → ∞ ∀ε >  <dig> where Λnε=1qn <dig> ∑k=1n∫|x|>εqnx2dvkxand qn2= ∑k=1nexk2and gx=2π∫0xe-t2/2dtif x ≥  <dig> and  <dig> if x <  <dig>  then we have supx|ps¯n<qnx-gx|→0where s¯n=max1≤k≤n ∑j=1kxjand vk = p
in order to apply the above theorem to get a simple formulaic approximation, we assume some random variables zi1m, each with the variance σ <dig> and sk= ∑i=1kzi.
then by applying the above theorem, we get the following uniform convergence of distribution to that of the one-sided standard normal as supx|pmaxk∈ <dig> …,msk≤mσx-gx|→0as m → ∞.
now we use the above results to get the probability estimates for our simple event terms {pij >xn} and {nij >xn}.
the following theorem provides us with the p-value's tail bound for lsa for any x∈ℝ.
theorem  <dig> for g, the one-sided normal distribution, defined above p|lsa|>x≤2n2-n-d-1n-d1-gxn/vαrx1y <dig>
proof.
by applying lemma  <dig> we have
ppij>xn=pmax <dig> pi- <dig> j-1+xi-1yj-1>xn=1-pmax <dig> pi- <dig> j-1+xi-1yj-1≤xn=1-pmax1≤k≤min{i- <dig> j-1}∑l=0kxlyl≤xn,
and by theorem  <dig>  replacing x with y, we have
supy|pmaxk∈ <dig> …,msk≤mσy-gy|→0asm→∞.
notice that ∑l=0kxlyl satisfies the definition of sk, so replacing sk, m , and σ with ∑l=0kxlyl, mini- <dig> j- <dig>  and varx1y <dig>  respectively,
supypmaxk∈ <dig> …,mini- <dig> j- <dig> ∑l=0kxlyl≤mini- <dig> j-1varx1y1y-gy→ <dig>
as n → ∞, and by change of variables to get our equation into the appropriate form
xn=mini- <dig> j-1varx1y1y⇒y=xn/mini- <dig> j-1varx1y1⇒supxpmaxkε <dig> …,mini- <dig> j- <dig> ∑l=0kxlyl≤xn-gxn/mini- <dig> j-1varx1y1→ <dig>
as n → ∞, thus
ppij>xn≅1-gxn/mini- <dig> j-1varx1y <dig>
it follows from boole's inequality and lemma  <dig> that
p|lsa|>x=p∪ijpij>xn∪∪ijnij>xn≤ ∑ijppij>xn+ ∑ijpnij>xn= <dig> ∑ij1-gxn/mini- <dig> j-1varx1y <dig>
finally, we have the following tail probability bound
p|lsa|>x≤ <dig> ∑ij1-gxn/mini- <dig> j-1varx1y1≤ <dig> ∑ij1-gxn/nvarx1y1=2n2-n-d-1n-d1-gxn/varx1y <dig>
standardizing with a mean of zero and a variance of one
=2n2-n-d-1n-d1-gxn.
□
note that this last result is asymptotic.
thus, n must be substantially large for this p-value bound to be relevant .
similar to the normal distribution as an approximation to student's t-distribution, this implementation of lsa requires at least  <dig> time points to promote confidence.
though this convergence can vary from dataset to dataset, the bound is conservative, and will not easily produce false positives if run on shorter time series.
RESULTS
to validate versatility and effectiveness of the derived upper bound , we applied the algorithm to four datasets, two sourced from biology, one from social networking, and a randomly generated control dataset.
these include the moving pictures of the human microbiome  <cit>  , the largest human microbial time series to date, a microarray hybridization dataset identifying cell cycle regulated genes in the yeast saccharomyces cerevisiae  <cit>  , and an online social media dataset of the volumes of the top  <dig> memetracker phrases and top  <dig> twitter hash tags over an eight month period from september  <dig> to august  <dig>  <cit> .
missing data values were interpolated by averaging the two nearest temporal data points, and all analysis was performed on a mac pro desktop computer running mac osx  <dig> .
<dig> with a  <dig> ×  <dig>  ghz quad-core intel xeon processors and  <dig> gb of  <dig> mhz ddr <dig> ram.
computational complexity
the algorithm calculates in o time, where m is the number of time series and n is the length of each time series.
to get an idea of how long calculations take, we fixed n =  <dig>  d =  <dig> and plotted log-calculation time against log-m .
it can be seen that lsa tests with p-values calculated by the permutation test are about  <dig>  times slower than calculating p-values formulaically.
compared to direct formulaic calculation, random number generation is slow, making a repetition of  <dig>  permutations for each time series pair a computationally intense operation .
the permutation test may be able to calculate statistically significant  pairs confidently, but applying a multiple test correction  will require exponentially more permutations to reach the same level of confidence for the entire dataset.
pairwise comparisons for big datasets are computationally infeasible to sufficiently estimate p-values with enough accuracy to protect against false positives.
in contrast, fastlsa directly calculates a conservative upper bound approximating the p-value, making permutation unnecessary and protecting against false positives.
moving pictures of the human microbiome
the mph time series dataset  <cit>  investigates temporal variations in the microbial community structure of two healthy human subjects, one male, one female.
samples were collected from three body parts, the gut , mouth, and skin  daily for  <dig> months  and six months  with taxonomy being determined by the amplified v <dig> region of the small subunit ribosomal rna  gene.
the male and female samples were concatenated together resulting in a profile of  <dig> taxa for  <dig> time points with missing values being interpolated by the average of the two nearest time points.
for a given time series, if more than 25% of time steps were zero it was removed from the analysis.
analysis took  <dig> minutes  without including output writing time which is variable.
significant  lsa results revealed clusters of local similarity that corresponded well when nodes were colored by sample source .
the low level of mixing between local clusters reflects the large differences in taxonomic profiles across the different body environments  <cit> .
microarray hybridization detection cell cycle-regulated genes in yeast saccharomyces cerevisae
in the cdc data set  <cit> , we focused on the cdc <dig> temperature sensitive mutant and the profile of  <dig> genes over  <dig> time steps, representing gene expression for approximately three cell cycles.
analysis took  <dig>  minutes  without including output writing time .
applying the asymptotic bound with the small number of time steps resulted in some rather large bounds .
however, lsa was capable of detecting lead-lag correlation despite the periodicity of the data, demonstrating its capacity to find long correlate pairs with a large number of covariate time series.
only  <dig> of the  <dig> gene nodes could be classified from  <cit>  to one of the five defined cell cycle phases  so only two clusters could be inferred upon with any confidence .
social media: top  <dig> twitter and memetracker phrases
the data from  <cit>  contains the volume of the top  <dig> twitter and memetracker phrase counts over  <dig> time steps from september  <dig> to august  <dig>  a spacing of approximately 2- <dig> days per observation.
analysis took approximately six seconds  without including write out time.
two major clusters of related times series nodes emerged.
however, attempts to label the series using existing metadata of general content or time granularity  did not elucidate its structure .
we conjecture that this difference is geographical  or socially structured, however, additional metadata on geolocation or social connectivity associations of the nodes would be needed to better elucidate network structure.
null hypothesis simulated data
finally, to identify throughput limits of fastlsa and to simulate a large iid dataset without time dependence, three data matrices were randomly generated:  one hundred thousand measurements across  <dig> time steps,  one million measurements across  <dig> time steps, and  one million measurements across  <dig> time steps.
data were generated by random sampling from a uniform distribution.
running fastlsa on  <dig> threads, the first dataset  took one hour  <dig> minutes, the second   <dig> days and  <dig> hours, and the third  had an eta of  <dig> days and  <dig> hours without including writeout time.
the asymptotic bound is conservative for shorter datasets  , the second data having  <dig> time points found zero false positives, despite having a bonferroni correction of αn2=10- <dig>  this is likely because the software's p-value is an upper bound to the real p-value, and so is the bonferroni correction.
an inspection of a uniform random graph  of  <dig>  random time series with  <dig> time steps did not generate any cliques, but only short  length chains of nodes, serving as a warning to those wanting to interpret relevant structure .
given appropriate thresholds on lsa values, cliques do not seem to occur randomly.
discussion
fastlsa uses a novel asymptotic upper bound algorithm for calculating the lsa p-value.
this is done without any normality assumption, extending implementation to untransformed data and data in violation of normality assumptions such as time series containing many zero entries.
moreover, fastlsa replaces a computationally intensive permutation test that was previously required to calculate significance of lsa statistics with a dramatic increase on the size of datasets that can be analyzed on a single desktop machine.
however, like all asymptotic bounds, a significant number of observations need to be obtained for their application.
from theoretical simulation, we estimate this to be greater than  <dig> time points for most datasets.
this is supported by our experience on the cdc and mph datasets having  <dig> and  <dig> time series, respectively.
despite this potential operating constraint, fastlsa expands the boundaries of lsa allowing time series analysis on datasets with millions of co-variate time series.
the algorithm is implemented as a computationally efficient software package, fastlsa, written in c and optimized for threading on multi-core computers using posix threads.
finally, we demonstrated the utility and versatility of fastlsa using real-world and simulated time series datasets from different fields of inquiry, visualizing the resulting clusters of local similarity using the cytoscape software.
lsa statistics have been demonstrated to capture relevant local similarity structure for a number of biological datasets  <cit> .
however, previous implementations were limited to relatively small datasets.
fastlsa improves the computational efficiency and statistical robustness of lsa, a necessary step in using the statistic to explore next generation time series datasets.
despite the current improvements, the structure captured by lsa is less than ideal and could be further improved.
given two vectors of time series, lsa reports the strongest statistic.
however, it is unclear where this significant time window occurs, or if there are multiple small windows with large lsa values that are not reported.
an inspection of time series traces in question is often required to visually check exactly how the two are similar.
another hazard is that lsa does not handle missing data effectively, and so a continuous version of the statistic would be desirable for exploratory experiments where sampling conditions could change to small degrees and analysis could be performed without imputation.
furthermore, lsa is asymmetric in nature, meaning that time reversal has the potential to produce differing lsa values.
we anticipate even better performance from the statistic if these issues were addressed, perhaps through a modified version of pcc that isolates optimal windows of similarity.
CONCLUSIONS
lsa is a local similarity statistic that has recently been used to capture relevant local structure in time series datasets, particularly within the biological community.
however, its use has been limited to smaller datasets due to an intensive permutation test used to calculate significance.
our derivation and direct calculation of an asymptotic upper bound using fastlsa replaces this onerous calculation without a normality assumption, enabling lsa on time series datasets containing millions of co-variate time series.
we demonstrate the utility and versatility of fastlsa by analyzing time series data from public health, microbial ecology, and social media and compare these results to the previous implementation of lsa, obtaining similar results with orders of magnitude increase in throughput.
project name: fastlsa
project home page: http://www.cmde.science.ubc.ca/hallam/fastlsa/
operating system: os x, linux, or windows
programming languages: c /c++
other requirements:  <dig> gb ram
license: gplv3
non-academic restrictions: none
list of abbreviations
lsa: local similarity analysis; pcc: pearson's correlation coefficient; pca: principal component analysis; mds: multidimensional scaling; dfa: discriminant fraction analysis; mph: moving pictures of the human microbiome; cdc: centre of disease control.
competing interests
the authors declare that they have no competing interests.
authors' contributions
wed derived the p-value upper bound and was the primary programmer of the lsa statistics.
nh assisted in the derivation of the upper bound, the analysis of the mph, cdc, and twitter datasets, and was the primary manuscript writer and editor.
kk validated upper bound result and assisted in the multi-threaded implementation of the software.
sh oversaw the research and managed the group.
declarations
the publication costs for this article were funded by genome british columbia and genome canada.
this article has been published as part of bmc genomics volume  <dig> supplement  <dig>  2013: selected articles from the eleventh asia pacific bioinformatics conference : genomics.
the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/14/s <dig>

 the exponential increase and ubiquitous use of computational technology has given rise to an era of "big data" that pushes the limits of conventional data analysis. techniques for analyzing big datasets often proceed by identifying patterns of co-occurrence or correlation through principal component analysis, multidimensional scaling, etc. however, many of these methods require significant data reduction or smoothing which makes them difficult to interpret. other methods such as multiple linear regression or pearson's correlation coefficient are easy to interpret as they operate on data in their native data space, without any kind of large data transformation or dimensionality reduction, but are limited in the structure that they can detect. though pcc is a classic and powerful technique for finding linear relationships between two variables, it is not designed for capturing lead-lag relationships seen in time series data. local similarity analysis extends correlation calculations to include the time variable, enabling identification of local correlates. furthermore, ruan et al. have presented a graphical network framework in which to visualize the structure of significant lsa correlations. unfortunately, the current implementation of lsa requires multiple runs on permuted data and a monte carlo statistical method known as a permutation test to evaluate a null distribution and obtain a p-value determining significance. each iteration of this procedure has a computational complexity of o, where p is the number of permutations, m is the number of covariate time series, and n is their length. due to the number of pair-wise calculations needed, extant lsa is computationally onerous when m is large, limiting its use to datasets where the number of observed variables at each time point is small. though there has been some improvement to its performance, assumptions of normality and implementation issues continue to stymie practical application of lsa on big datasets. here we describe a novel asymptotic upper bound on the calculation of the lsa statistic's p-value, resulting in an exponentially converging calculation to bound and check for significance of computed lsa statistics without a computationally intensive permutation test. this bound does not require a rank-vector normal transformation, promoting application to any distribution that has finite variance. as a result, this implementation of lsa can navigate big datasets with millions of co-variate time series. we demonstrate this using time-series datasets from public health, microbial ecology, and social media. the implemented algorithm, named fastlsa, is written in c and optimized for threading on multi-core computers. interpreting the lsa statistic lsa concerns itself with pairs of time series data. the lsa statistic can be interpreted in a manner similar to pcc when no lag window exists between two time series. however, lsa is also capable of capturing localized correlation that is staggered or lagged. a large positive or negative lsa value indicates a correspondingly strong pcc correlation or a correlation at a time displacement within the lag window. note that if both positive and negative correlations exist between two time series, lsa will only report the strongest of the two. lsa is advantageous on large datasets containing many time series. results can be visualized as a graphical network where nodes represent the individual time series and the edges represent their lsa correlation statistic. when displayed using a force-directed layout in cytoscape, closely related time series cluster together, visually isolating clusters of local similarity. metadata related to experimental or environmental conditions can then be applied to the nodes, shedding insight into hierarchical network structure. implementation description of the lsa algorithm in this section we reproduce the algorithm from to compute lsa statistics and their corresponding p-values between pairs of time series in a dataset. we assume as input a set of time series vectors of equal length. let us denote the number of time series by m and their length as n. let us denote the time series dataset as x where xij denotes the jth element of the ith time series, with i =..., m and j =..., n, and assume that the xij are real numbers. we also assume that there are no missing values in the dataset x, and realize that practical use will require interpolation or filtering. in figure we present the algorithm for computing the lsa statistic for a pair of time series, x=xi1n and y=yi1n, where the length of the time series is assumed to be equally spaced in time. we have modified the presentation of the lsa algorithm by to highlight our analysis and derivation of a bound on the tail distribution of the lsa statistic. specifically we calculate the lsa statistic for a pair of time series, x and y. two-dimensional arrays pi, j and ni, j are used to store the positive and negative partial sums of the pairwise product of time series values. we also assume a suitable bound on the maximum time lag considered while computing the lsa statistic, denoted by d. the algorithm first initializes the arrays pj, nj, p i, and n i for all i, j =..., n, with a maximum absolute difference of d. next it considers the time series pairs for each possible lag, up to a maximum of d, and then computes the progressive sum of the pair-wise products of the time-series values from the low to high index of the arrays. during the computation, the progression of the partial sum is reset to if the sum is below after partial sums have been computed, the values of n ^ and p ^ are calculated by taking the maximum of the corresponding values of the arrays n and p. finally, the lsa statistic is estimated as signp ^-n ^maxp ^,n ^n. calculating the upper bound in this section we derive the asymptotic upper bound on the p-value for the cumulative probability distribution of the lsa statistic without the need of a normality assumption. our derivation is based on distributional results of the maximum cumulative sum of independent random variables known in the literature from probability theory. we begin by stating our assumptions about the dataset, isolate target calculations from the lsa algorithm, and from our referenced mathematical results, derive and prove important lemmas. these lemmas will serve as the building blocks as we logically construct a theorem which will form the basis of our lsa p-value upper bound. we begin by making certain assumptions about the probability model used to derive the bounds. first, each pi, j or ni, j is considered individually. we assume that the time series values xi, yj for i, j =..., n are independent of one another. this assumption can be made when weak dependence exists because it is near the truth and effective, much like the naive bayes assumption. this assumption is also enabling, as it allows us to invoke the distributions of partial sums of independent random variables and continue in a mathematically straightforward way. further, we assume independence between each time time series as a null hypothesis, and as it is subject to rejection upon obtaining a statistically significant lsa value. consider lines and of the lsa algorithm, pi+k+ j+k+ ← max{ pi+k, j+k + xi+k * yj+k} and p ^←maxi,j:|i-j|≤dpi,j. for any pair of i and j let us define the sequence random variables as zk = xi+kyj+k for k =..., min{n - i, n - j} - and the sequence of random variables ζk = z +... + zk for k =..., min{n - i, n - j} - supposing ζ = using the above ζk's, we define random variables ηk* as ηk*=maxζ ζ ⋯,ζk for the same values of k =..., min{n - i, n - j} - we also define the set of random variables η η..., ηk by the recurrence formula ηk+ = max{ ηk + zk+1}. note that the random variables pi+k, j+k and the ηk have the same distribution. it is shown in that the random variables ηk* and ηk also have the same distribution. as a result, now we can analyze the cumulative distribution of pi+k, j+k as a distribution for ηk*, and use the results by nevzorov and petrov on pi+k, j+k to derive tail probability bounds. we also assume that the random variables zk have the first two moments, although such assumptions are not required for the results of, we use them to derive simpler bounds. we now consider a few useful lemmas that we will use to construct our p-value upper bound. the first step is to simplify the tail event into simpler terms. the following lemma expresses the tail event for lsa {|lsa| >x} and any x∈ℝ in terms of the tail events of {pi, j >x} and {ni, j >x}, the positive and negative lsa calculations for the same x, the bound on our test statistic. lemma for any x∈ℝwe have {|lsa| >x} = { ∪ }. proof. the result is clear from the following: |lsa|>x=maxp ^,n ^>xn=p ^≤xn∩n ^≤xnc=maxijpij≤xn∩maxijnij≤xnc=∩ijpij≤xn∩∩ijnij≤xnc=∪ijpij>xn∪∪ijnij>xn □ in the lsa algorithm, we have maximums pij = max{ pi- j- + xi-1yj-1} and nij = max{ pi- j- - xi- j-1}, which complicates their theoretical analysis. fortunately, equivalence have been demonstrated in the literature, and we restate these in the following lemma for clarity: the similarity of the distributions of ηk and ηk*, for k =..., min{n - i, n - j} - this will help us derive the bounds for the events {pij >xn} and {nij >xn}, the simpler terms we derived in the previous lemma. lemma let zi be mutually independent random variables and let us denote by sk= ∑i=1kziwhere s = and qk+ = max{ qk + zk} with q = then p = p for x∈ℝ. in order to get a simple formula for the bound on the cumulative tail probabilities for pi, j and ni, j we reproduce below the results on partial sums of random variables due to nevzorov and petrov. for our sequence of independent and identically distributed random variables under consideration {xn} it follows that lindeberg's condition holds. a property showing the variance of a distribution stabilizes as more variables are added, pinning the tails of it down. thinking about this in terms of time series, as a series gets larger, the upper bound of the distribution becomes more defined and calculable. now to build theorems upon which we will derive a formulaic p-value bound. theorem if the random variables {xn} have zero expectation and finite variances and if lindeberg's condition holds: Λn → as n → ∞ ∀ε > where Λnε=1qn ∑k=1n∫|x|>εqnx2dvkxand qn2= ∑k=1nexk2and gx=2π∫0xe-t2/2dtif x ≥ and if x < then we have supx|ps¯n<qnx-gx|→0where s¯n=max1≤k≤n ∑j=1kxjand vk = p in order to apply the above theorem to get a simple formulaic approximation, we assume some random variables zi1m, each with the variance σ and sk= ∑i=1kzi. then by applying the above theorem, we get the following uniform convergence of distribution to that of the one-sided standard normal as supx|pmaxk∈ …,msk≤mσx-gx|→0as m → ∞. now we use the above results to get the probability estimates for our simple event terms {pij >xn} and {nij >xn}. the following theorem provides us with the p-value's tail bound for lsa for any x∈ℝ. theorem for g, the one-sided normal distribution, defined above p|lsa|>x≤2n2-n-d-1n-d1-gxn/vαrx1y proof. by applying lemma we have ppij>xn=pmax pi- j-1+xi-1yj-1>xn=1-pmax pi- j-1+xi-1yj-1≤xn=1-pmax1≤k≤min{i- j-1}∑l=0kxlyl≤xn, and by theorem replacing x with y, we have supy|pmaxk∈ …,msk≤mσy-gy|→0asm→∞. notice that ∑l=0kxlyl satisfies the definition of sk, so replacing sk, m, and σ with ∑l=0kxlyl, mini- j- and varx1y respectively, supypmaxk∈ …,mini- j- ∑l=0kxlyl≤mini- j-1varx1y1y-gy→ as n → ∞, and by change of variables to get our equation into the appropriate form xn=mini- j-1varx1y1y⇒y=xn/mini- j-1varx1y1⇒supxpmaxkε …,mini- j- ∑l=0kxlyl≤xn-gxn/mini- j-1varx1y1→ as n → ∞, thus ppij>xn≅1-gxn/mini- j-1varx1y it follows from boole's inequality and lemma that p|lsa|>x=p∪ijpij>xn∪∪ijnij>xn≤ ∑ijppij>xn+ ∑ijpnij>xn= ∑ij1-gxn/mini- j-1varx1y finally, we have the following tail probability bound p|lsa|>x≤ ∑ij1-gxn/mini- j-1varx1y1≤ ∑ij1-gxn/nvarx1y1=2n2-n-d-1n-d1-gxn/varx1y standardizing with a mean of zero and a variance of one =2n2-n-d-1n-d1-gxn. □ note that this last result is asymptotic. thus, n must be substantially large for this p-value bound to be relevant. similar to the normal distribution as an approximation to student's t-distribution, this implementation of lsa requires at least time points to promote confidence. though this convergence can vary from dataset to dataset, the bound is conservative, and will not easily produce false positives if run on shorter time series.  to validate versatility and effectiveness of the derived upper bound, we applied the algorithm to four datasets, two sourced from biology, one from social networking, and a randomly generated control dataset. these include the moving pictures of the human microbiome, the largest human microbial time series to date, a microarray hybridization dataset identifying cell cycle regulated genes in the yeast saccharomyces cerevisiae, and an online social media dataset of the volumes of the top memetracker phrases and top twitter hash tags over an eight month period from september to august. missing data values were interpolated by averaging the two nearest temporal data points, and all analysis was performed on a mac pro desktop computer running mac osx. with a × ghz quad-core intel xeon processors and gb of mhz ddr ram. computational complexity the algorithm calculates in o time, where m is the number of time series and n is the length of each time series. to get an idea of how long calculations take, we fixed n = d = and plotted log-calculation time against log-m. it can be seen that lsa tests with p-values calculated by the permutation test are about times slower than calculating p-values formulaically. compared to direct formulaic calculation, random number generation is slow, making a repetition of permutations for each time series pair a computationally intense operation. the permutation test may be able to calculate statistically significant pairs confidently, but applying a multiple test correction will require exponentially more permutations to reach the same level of confidence for the entire dataset. pairwise comparisons for big datasets are computationally infeasible to sufficiently estimate p-values with enough accuracy to protect against false positives. in contrast, fastlsa directly calculates a conservative upper bound approximating the p-value, making permutation unnecessary and protecting against false positives. moving pictures of the human microbiome the mph time series dataset investigates temporal variations in the microbial community structure of two healthy human subjects, one male, one female. samples were collected from three body parts, the gut, mouth, and skin daily for months and six months with taxonomy being determined by the amplified v region of the small subunit ribosomal rna gene. the male and female samples were concatenated together resulting in a profile of taxa for time points with missing values being interpolated by the average of the two nearest time points. for a given time series, if more than 25% of time steps were zero it was removed from the analysis. analysis took minutes without including output writing time which is variable. significant lsa results revealed clusters of local similarity that corresponded well when nodes were colored by sample source. the low level of mixing between local clusters reflects the large differences in taxonomic profiles across the different body environments. microarray hybridization detection cell cycle-regulated genes in yeast saccharomyces cerevisae in the cdc data set, we focused on the cdc temperature sensitive mutant and the profile of genes over time steps, representing gene expression for approximately three cell cycles. analysis took minutes without including output writing time. applying the asymptotic bound with the small number of time steps resulted in some rather large bounds. however, lsa was capable of detecting lead-lag correlation despite the periodicity of the data, demonstrating its capacity to find long correlate pairs with a large number of covariate time series. only of the gene nodes could be classified from to one of the five defined cell cycle phases so only two clusters could be inferred upon with any confidence. social media: top twitter and memetracker phrases the data from contains the volume of the top twitter and memetracker phrase counts over time steps from september to august a spacing of approximately 2- days per observation. analysis took approximately six seconds without including write out time. two major clusters of related times series nodes emerged. however, attempts to label the series using existing metadata of general content or time granularity did not elucidate its structure. we conjecture that this difference is geographical or socially structured, however, additional metadata on geolocation or social connectivity associations of the nodes would be needed to better elucidate network structure. null hypothesis simulated data finally, to identify throughput limits of fastlsa and to simulate a large iid dataset without time dependence, three data matrices were randomly generated: one hundred thousand measurements across time steps, one million measurements across time steps, and one million measurements across time steps. data were generated by random sampling from a uniform distribution. running fastlsa on threads, the first dataset took one hour minutes, the second days and hours, and the third had an eta of days and hours without including writeout time. the asymptotic bound is conservative for shorter datasets, the second data having time points found zero false positives, despite having a bonferroni correction of αn2=10- this is likely because the software's p-value is an upper bound to the real p-value, and so is the bonferroni correction. an inspection of a uniform random graph of random time series with time steps did not generate any cliques, but only short length chains of nodes, serving as a warning to those wanting to interpret relevant structure. given appropriate thresholds on lsa values, cliques do not seem to occur randomly. discussion fastlsa uses a novel asymptotic upper bound algorithm for calculating the lsa p-value. this is done without any normality assumption, extending implementation to untransformed data and data in violation of normality assumptions such as time series containing many zero entries. moreover, fastlsa replaces a computationally intensive permutation test that was previously required to calculate significance of lsa statistics with a dramatic increase on the size of datasets that can be analyzed on a single desktop machine. however, like all asymptotic bounds, a significant number of observations need to be obtained for their application. from theoretical simulation, we estimate this to be greater than time points for most datasets. this is supported by our experience on the cdc and mph datasets having and time series, respectively. despite this potential operating constraint, fastlsa expands the boundaries of lsa allowing time series analysis on datasets with millions of co-variate time series. the algorithm is implemented as a computationally efficient software package, fastlsa, written in c and optimized for threading on multi-core computers using posix threads. finally, we demonstrated the utility and versatility of fastlsa using real-world and simulated time series datasets from different fields of inquiry, visualizing the resulting clusters of local similarity using the cytoscape software. lsa statistics have been demonstrated to capture relevant local similarity structure for a number of biological datasets. however, previous implementations were limited to relatively small datasets. fastlsa improves the computational efficiency and statistical robustness of lsa, a necessary step in using the statistic to explore next generation time series datasets. despite the current improvements, the structure captured by lsa is less than ideal and could be further improved. given two vectors of time series, lsa reports the strongest statistic. however, it is unclear where this significant time window occurs, or if there are multiple small windows with large lsa values that are not reported. an inspection of time series traces in question is often required to visually check exactly how the two are similar. another hazard is that lsa does not handle missing data effectively, and so a continuous version of the statistic would be desirable for exploratory experiments where sampling conditions could change to small degrees and analysis could be performed without imputation. furthermore, lsa is asymmetric in nature, meaning that time reversal has the potential to produce differing lsa values. we anticipate even better performance from the statistic if these issues were addressed, perhaps through a modified version of pcc that isolates optimal windows of similarity.  lsa is a local similarity statistic that has recently been used to capture relevant local structure in time series datasets, particularly within the biological community. however, its use has been limited to smaller datasets due to an intensive permutation test used to calculate significance. our derivation and direct calculation of an asymptotic upper bound using fastlsa replaces this onerous calculation without a normality assumption, enabling lsa on time series datasets containing millions of co-variate time series. we demonstrate the utility and versatility of fastlsa by analyzing time series data from public health, microbial ecology, and social media and compare these results to the previous implementation of lsa, obtaining similar results with orders of magnitude increase in throughput. project name: fastlsa project home page: http://www.cmde.science.ubc.ca/hallam/fastlsa/ operating system: os x, linux, or windows programming languages: c /c++ other requirements: gb ram license: gplv3 non-academic restrictions: none list of abbreviations lsa: local similarity analysis; pcc: pearson's correlation coefficient; pca: principal component analysis; mds: multidimensional scaling; dfa: discriminant fraction analysis; mph: moving pictures of the human microbiome; cdc: centre of disease control. competing interests the authors declare that they have no competing interests. authors' contributions wed derived the p-value upper bound and was the primary programmer of the lsa statistics. nh assisted in the derivation of the upper bound, the analysis of the mph, cdc, and twitter datasets, and was the primary manuscript writer and editor. kk validated upper bound result and assisted in the multi-threaded implementation of the software. sh oversaw the research and managed the group. declarations the publication costs for this article were funded by genome british columbia and genome canada. this article has been published as part of bmc genomics volume supplement 2013: selected articles from the eleventh asia pacific bioinformatics conference: genomics. the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/14/s 
we computationally compare our approach to previous implementations of lsa, demonstrate broad applicability by analyzing time series data from public health, microbial ecology, and social media, and visualize resulting networks using the cytoscape software. pairwise comparison of time series data for both local and time-lagged relationships is a computationally challenging problem relevant to many fields of inquiry. to improve the performance of lsa on big datasets, an asymptotic upper bound on the p-value calculation was derived without the assumption of normality. the local similarity analysis statistic identifies the existence of local and lagged relationships, but determining significance through a p-value has been algorithmically cumbersome due to an intensive permutation test, shuffling rows and columns and repeatedly calculating the statistic. this change in the bound calculation markedly improved computational speed from o to o, where p is the number of permutations in a permutation test, m is the number of time series, and n is the length of each time series. the fastlsa software package expands the boundaries of lsa allowing analysis on datasets with millions of co-varying time series. the bounding process is implemented as a computationally efficient software package, fastlsa, written in c and optimized for threading on multi-core computers, improving its practical computation time.
 pairwise comparison of time series data for both local and time-lagged relationships is a computationally challenging problem relevant to many fields of inquiry. the local similarity analysis statistic identifies the existence of local and lagged relationships, but determining significance through a p-value has been algorithmically cumbersome due to an intensive permutation test, shuffling rows and columns and repeatedly calculating the statistic. furthermore, this p-value is calculated with the assumption of normality -- a statistical luxury dissociated from most real world datasets.  to improve the performance of lsa on big datasets, an asymptotic upper bound on the p-value calculation was derived without the assumption of normality. this change in the bound calculation markedly improved computational speed from o to o, where p is the number of permutations in a permutation test, m is the number of time series, and n is the length of each time series. the bounding process is implemented as a computationally efficient software package, fastlsa, written in c and optimized for threading on multi-core computers, improving its practical computation time. we computationally compare our approach to previous implementations of lsa, demonstrate broad applicability by analyzing time series data from public health, microbial ecology, and social media, and visualize resulting networks using the cytoscape software.  the fastlsa software package expands the boundaries of lsa allowing analysis on datasets with millions of co-varying time series. mapping metadata onto force-directed graphs derived from fastlsa allows investigators to view correlated cliques and explore previously unrecognized network relationships. the software is freely available for download at: http://www.cmde.science.ubc.ca/hallam/fastlsa/. 21- january the eleventh asia pacific bioinformatics conference vancouver, canada 