# Online algorithm for mean absolute deviation and large data set

I have a little problem that is making me freaking out. I have to write procedure for an online acquisition process of a multivariate time series. At every time interval (for example 1 second), I get a new sample, which is basically a floating point vector of size N. The operation I need to do is a little bit tricky:

1. For each new sample, I compute the percentuals for that sample (by normalizing the vector so that the elements will sum to 1).

2. I calculate the average percentuals vector in the same way, but using the past values.

3. For each past value, I compute the absolute deviation of the percentuals vector related to that sample with the global average percentuals vector computed at step 2. This way, the absolute deviation is always a number between 0 (when the vector is equal to the average vector) and 2 (when it is totaly different).

4. Using the average of the deviations for all the previous samples, I compute the mean absolute deviation, which is again a number between 0 and 2.

5. I use the mean absolute deviation to detect if a new sample is compatible with the other samples (by comparing its absolute deviation with the mean absolute deviation of the whole set computed at step 4).

Since every time a new sample is collected the global average changes (and so the mean absolute deviation changes as well), is there a way to compute this value without scanning the whole data set multiple times? (one time for computing the global average percentuals, and one time for collecting the absolute deviations). Ok, I know it's absolutely easy to calculate the global averages without scanning the whole set, since I just have to use a temporary vector for storing the sum of each dimension, but what about the mean absolute deviation? Its calculation includes the abs() operator, so I need to access to all the past data!

-

If you can accept some inaccuracy, this problem can be solved easily by binning counts. That is, pick some largeish number $M$ (say $M = 1000$), then initialize some integer bins $B_{i,j}$ for $i = 1\ldots M$ and $j = 1\ldots N$, where $N$ is the vector size, as zero. Then when you see the $k$th observation of a percentual vector, increment $B_{i,j}$ if the $j$th element of this vector is between $(i-1)/M$ and $i/M$, looping over the $N$ elements of the vector. (I am assuming your input vectors are non-negative, so that when you compute your 'percentuals', the vectors are in the range $[0,1]$. )

At any point in time, you can estimate the mean vector from the bins, and the mean absolute deviation. After observing $K$ such vectors, the $j$th element of the mean is estimated by $$\bar{X}_j = \frac{1}{K} \sum_i \frac{i - 1/2}{M} B_{i,j},$$ and the $j$th element of the mean absolute deviation is estimated by $$\frac{1}{K} \sum_i | \bar{X_j} - \frac{i - 1/2}{M} | B_{i,j}$$

edit: this is a specific case of a more general approach where you are building an empirical density estimate. This could be done with polynomials, splines, etc, but the binning approach is the easiest to describe and implement.

-
Wow, really interesting approach. I didn't know about that, and I will keep it in mind. Unfortunately, in this case it will not work, since I have really restrictive requirements from the point of view of memory usage, so M should be really small, and I guess there would be too much precision loss. –  gianluca Oct 7 '10 at 5:26
@gianluca: it sounds like you have 1. a lot of data, 2. limited memory resources, 3. high precision requirements. I can see why this problem is freaking you out! Perhaps, as mentioned by @kwak, you can compute some other measure of spread: MAD, IQR, standard deviation. All of those have approaches which might work for your problem. –  shabbychef Oct 7 '10 at 16:22
gianluca:> Give us more quantitative idea about the size of memory, arrays and accuracy you want. It may well be that your question will be best answered @ stackoverflow though. –  user603 Oct 7 '10 at 18:28

I've used the following approach in the past to calculate absolution deviation moderately efficiently (note, this a programmers approach, not a statisticians, so indubitably there may be clever tricks like shabbychef's that might be more efficient).

WARNING: This is not an online algorithm. It requires O(n) memory. Furthermore, it has a worst case performance of O(n), for datasets like [1, -2, 4, -8, 16, -32, ...] (i.e. the same as the full recalculation). [1]

However, because it still performs well in many use cases it might be worth posting here. For example, in order to calculate the absolute deviance of 10000 random numbers between -100 and 100 as each item arrives, my algorithm takes less than one second, while the full recalculation takes over 17 seconds (on my machine, will vary per machine and according to input data). You need to maintain the entire vector in memory however, which may be a constraint for some uses. The outline of the algorithm is as follows:

1. Instead of having a single vector to store past measurements, use three sorted priority queues (something like a min/max heap). These three lists partition the input into three: items greater than the mean, items less than the mean and items equal to the mean.
2. (Almost) every time you add an item the mean changes, so we need to repartition. The crucial thing is the sorted nature of the partitions which means that instead of scanning every item in the list to repartion, we only need to read those items we are moving. While in the worst case this will still require O(n) move operations, for many use-cases this is not so.
3. Using some clever bookkeeping, we can make sure that the deviance is correctly calculated at all times, when repartitioning and when adding new items.

Some sample code, in python, is below. Note that it only allows items to be added to the list, not removed. This could easily be added, but at the time I wrote this I had no need for it. Rather than implement the priority queues myself, I have used the sortedlist from Daniel Stutzbach's excellent blist package, which use B+Trees internally.

Consider this code licensed under the MIT license. It has not been significantly optimised or polished, but has worked for me in the past. New versions will be available here. Let me know if you have any questions, or find any bugs.

from blist import sortedlist
import operator

class deviance_list:
def __init__(self):
self.mean =  0.0
self._old_mean = 0.0
self._sum =  0L
self._n =  0  #n items
# items greater than the mean
self._toplist =  sortedlist()
# items less than the mean
self._bottomlist = sortedlist(key = operator.neg)
# Since all items in the "eq list" have the same value (self.mean) we don't need
# to maintain an eq list, only a count
self._eqlistlen = 0

self._top_deviance =  0
self._bottom_deviance =  0

@property
def absolute_deviance(self):
return self._top_deviance + self._bottom_deviance

def append(self,  n):
# Update summary stats
self._sum += n
self._n +=  1
self._old_mean =  self.mean
self.mean =  self._sum /  float(self._n)

# Move existing things around
going_up = self.mean > self._old_mean
self._rebalance(going_up)

# Add new item to appropriate list
if n >  self.mean:
self._top_deviance +=  n -  self.mean
elif n == self.mean:
self._eqlistlen += 1
else:
self._bottom_deviance += self.mean -  n

def _move_eqs(self,  going_up):
if going_up:
self._bottomlist.update([self._old_mean] *  self._eqlistlen)
self._bottom_deviance += (self.mean - self._old_mean) * self._eqlistlen
self._eqlistlen = 0
else:
self._toplist.update([self._old_mean] *  self._eqlistlen)
self._top_deviance += (self._old_mean - self.mean) * self._eqlistlen
self._eqlistlen = 0

def _rebalance(self, going_up):
move_count,  eq_move_count = 0, 0
if going_up:
# increase the bottom deviance of the items already in the bottomlist
if self.mean !=  self._old_mean:
self._bottom_deviance += len(self._bottomlist) *  (self.mean -  self._old_mean)
self._move_eqs(going_up)

# transfer items from top to bottom (or eq) list, and change the deviances
for n in iter(self._toplist):
if n < self.mean:
self._top_deviance -= n -  self._old_mean
self._bottom_deviance += (self.mean -  n)
# we increment movecount and move them after the list
# has finished iterating so we don't modify the list during iteration
move_count +=  1
elif n == self.mean:
self._top_deviance -= n -  self._old_mean
self._eqlistlen += 1
eq_move_count +=  1
else:
break
for _ in xrange(0,  move_count):
for _ in xrange(0,  eq_move_count):
self._toplist.pop(0)

# decrease the top deviance of the items remain in the toplist
self._top_deviance -= len(self._toplist) *  (self.mean -  self._old_mean)
else:
if self.mean !=  self._old_mean:
self._top_deviance += len(self._toplist) *  (self._old_mean -  self.mean)
self._move_eqs(going_up)
for n in iter(self._bottomlist):
if n > self.mean:
self._bottom_deviance -= self._old_mean -  n
self._top_deviance += n -  self.mean
move_count += 1
elif n == self.mean:
self._bottom_deviance -= self._old_mean -  n
self._eqlistlen += 1
eq_move_count +=  1
else:
break
for _ in xrange(0,  move_count):
for _ in xrange(0,  eq_move_count):
self._bottomlist.pop(0)

# decrease the bottom deviance of the items remain in the bottomlist
self._bottom_deviance -= len(self._bottomlist) *  (self._old_mean -  self.mean)

if __name__ ==  "__main__":
import random
dv =  deviance_list()
# Test against some random data,  and calculate result manually (nb. slowly) to ensure correctness
rands = [random.randint(-100,  100) for _ in range(0,  1000)]
ns = []
for n in rands:
dv.append(n)
ns.append(n)
(n, dv.mean,  dv._old_mean,  dv.absolute_deviance / dv.mean))
assert sum(ns) == dv._sum,  "Sums not equal!"
assert len(ns) == dv._n,  "Counts not equal!"
m = sum(ns) / float(len(ns))
assert m == dv.mean,  "Means not equal!"
real_abs_dev = sum([abs(m - x) for x in ns])
# Due to floating point imprecision, we check if the difference between the
# two ways of calculating the asb. dev. is small rather than checking equality
assert abs(real_abs_dev - dv.absolute_deviance) < 0.01, (
"Absolute deviances not equal. Real:%.2f,  calc:%.2f" %  (real_abs_dev,  dv.absolute_deviance))


[1] If symptoms persist, see your doctor.

-
I'm missing something: if you have to "maintain the entire vector in memory," how does this qualify as an "online" algorithm?? –  whuber Nov 6 '10 at 15:03
@whuber No, not missing something, I guess it isn't an online algorithm. It requires O(n) memory, and in the worst case takes O(n) time for each item added. In normally distributed data (and probably other distributions) it works quite efficiently though. –  fmark Nov 9 '10 at 5:29

There is also a parametric approach. Ignoring the vector nature of your data, and looking only at the marginals, it suffices to solve the problem: find an online algorithm to compute the mean absolute deviation of scalar $X$. If (and this is the big 'if' here) you thought that $X$ followed some probability distribution with unknown parameters, you could estimate the parameters using some online algorithm, then compute the mean absolute deviation based on that parametrized distribution. For example, if you thought that $X$ was (approximately) normally distributed, you could estimate its standard deviation, as $s$, and the mean absolute deviation would be estimated by $s \sqrt{2 / \pi}$ (see Half Normal Distribution).

-
That's an interesting idea. You could supplement it, perhaps, with online detection of outliers and use those to modify the estimate as you go along. –  whuber Nov 6 '10 at 15:05
You could probably use Welford's method for calculating the standard deviation online that I documented in my second answer. –  fmark Nov 9 '10 at 5:38
One should however note that this way one might lose robustness of estimators such as explicit MAD, which are sometimes driving to its choice against simpler alternatives. –  Quartz Sep 13 at 15:47

MAD(x) is just two concurrent median computation, each of which can be made online through the binmedian algorithm.

You can find the associated paper as well as C and FORTRAN code online here.

(this is just the use of a clever trick on top of Shabbychef's clever trick, to save on memory).

There exist a host of older multi-pass methods for computing quantiles. A popular approach is to maintain/update a deterministically sized reservoir of observations randomly selected from the stream and recursively compute quantiles (see this review) on this reservoir. This (and related) approach are superseded by the one proposed above.

-
Could you please detail or reference the relation between MAD and the two medians? –  Quartz Sep 13 at 15:47
it's really the formula of the MAD: $\text{med}_{i=1}^n|x_i-\text{med}_{i=1}^n|$ (hence two medians) –  user603 Sep 14 at 16:02
Hehm, I actually meant if you can explain how is this relation allowing for the two medians to be concurrent; those seem dependent to me, since the inputs to the outer median may all change at each added sample to the inner calculation. How would you perform them in parallel? –  Quartz Sep 15 at 13:31
I have to go back to the binmedian paper for details...but given a computed value of the median ($med_{i=1}^n x_i$) and a new value of $x_{n+1}$ the algorithm could compute $med_{i=1}^{n+1}x_i$ much faster than $\mathcal{O}(n)$ by identifying the bin to which $x_{n+1}$ belongs. I don't see how this insight could not be generalized to the outer median in the mad computation. –  user603 Sep 15 at 13:55

The following provides an inaccurate approximation, although the inaccuracy will depend on the distribution of the input data. It is an online algorithm, but only approximates the absolute deviance. It is based on a well known algorithm for calculating variance online, described by Welford in the 1960s. His algorithm, translated into R, looks like:

M2 <- 0
mean <- 0
n <- 0

var.online <- function(x){
n <<- n + 1
diff <- x - mean
mean <<- mean + diff / n
M2 <<- M2 + diff * (x - mean)
variance <- M2 / (n - 1)
return(variance)
}


It performs very similarly to R's builtin variance function:

set.seed(2099)
n.testitems <- 1000
n.tests <- 100
differences <- rep(NA, n.tests)
for (i in 1:n.tests){
# Reset counters
M2 <- 0
mean <- 0
n <- 0

xs <- rnorm(n.testitems)
for (j in 1:n.testitems){
v <- var.online(xs[j])
}

differences[i] <- abs(v - var(xs))

}
summary(differences)
Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
0.000e+00 2.220e-16 4.996e-16 6.595e-16 9.992e-16 1.887e-15


Modifying the algorithm to calculate absolute deviation simply involves an additional sqrt call. However, the sqrt introduces inaccuracies that are reflected in the result:

absolute.deviance.online <- function(x){
n <<- n + 1
diff <- x - mean
mean <<- mean + diff / n
a.dev <<- a.dev + sqrt(diff * (x - mean))
return(a.dev)
}


The errors, calculated as above, are much greater than for the variance calculation:

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
0.005126 0.364600 0.808000 0.958800 1.360000 3.312000


However, depending on your use case, this magnitude of error might be acceptable.

-
This does not give the exact answer, for the following reason: $\sqrt{\sum_i x_i} \ne \sum_i \sqrt{x_i}$. You are computing the former, while the OP wants the latter. –  shabbychef Oct 12 '10 at 19:45
I agree that the method is inexact. However, I disagree with your diagnosis of the inexactness. Welford's method for calculating variance, which does not even contain a sqrt, has a similar error. However, as n gets large, the error/n gets vanishingly small, suprisingly quickly. –  fmark Oct 12 '10 at 23:06
Welford's method has no sqrt because it is computing the variance, not the standard deviation. By taking the sqrt, it seems like you are estimating the standard deviation, not the mean absolute deviation. am I missing something? –  shabbychef Oct 12 '10 at 23:33
@shabbychef Each iteration of Welfords is calculating the contribution of the new datapoint to the absolute deviation, squared. So I take the square root of each contribution squared to get back to the absolute deviance. You might note, for example, that I take the square root of the delta before I add it to the deviance sum, rather than afterward as in the case of the standard deviation. –  fmark Oct 13 '10 at 5:07
I see the problem; Welfords obscures the problem with this method: the online estimate of the mean is being used instead of the final estimate of the mean. While Welford's method is exact (up to roundoff) for variance, this method is not. The problem is not due to the sqrt imprecision. It is because it uses the running mean estimate. To see when this will break, try xs <- sort(rnorm(n.testitems)) When I try this with your code (after fixing it to return a.dev / n), I get relative errors on the order of 9%-16%. So this method is not permutation invariant, which could cause havoc... –  shabbychef Oct 13 '10 at 17:05