Running Standard Deviation in MySQL

MySQL's built-in STD() (an alias for STDDEV_POP()) computes the standard deviation of a column in a single aggregate pass. That works fine as long as the raw rows are still around, but it stops scaling the moment you want a live, per-group figure that updates as new rows arrive: recomputing STD() from scratch after every insert costs \(O(n)\) per group. If you don't need to keep the raw rows, it pays off to instead maintain a small per-group summary that a single new row can update in \(O(1)\).

Setting Up a Baseline

A test table with some random values, grouped into five buckets, gives something to compare against:

CREATE TABLE t1(
  group_value INT UNSIGNED NOT NULL,
  value DOUBLE NOT NULL
);

INSERT INTO t1 (group_value, value)
SELECT FLOOR(RAND() * 5), RAND() * 1000
FROM information_schema.tables
LIMIT 50;

The ordinary aggregate query gives the baseline every later approach has to reproduce:

SELECT group_value, COUNT(*), AVG(value), STD(value)
FROM t1
GROUP BY group_value;

+-------------+----------+--------------------+--------------------+
| group_value | COUNT(*) | AVG(value)         | STD(value)         |
+-------------+----------+--------------------+--------------------+
|           0 |       13 |  558.1130023798883 |    307.22234586301 |
|           1 |       10 |  560.2472776176849 |  287.3629941826582 |
|           2 |       11 | 426.06415503799786 | 219.53149074497404 |
|           3 |        7 |  486.4831761597233 | 315.31028664538627 |
|           4 |        9 |  540.8957763955889 |  221.4856395901883 |
+-------------+----------+--------------------+--------------------+
5 rows in set (0.00 sec)

STD() is the population standard deviation, i.e. it divides by \(n\), not \(n-1\). Every formula derived below targets that same population form, so it stays comparable to this baseline.

Sum and Sum of Squares

The population variance of a sample \(x_1,\dots,x_n\) with mean \(\bar x\) is

\[ \sigma^2 = \frac1n\sum_{i=1}^n (x_i-\bar x)^2. \]

Expanding the square and using \(\bar x=\frac1n\sum x_i\) removes every reference to the individual deviations:

\[ \begin{array}{rl} \sigma^2 &= \dfrac1n\sum x_i^2 - \dfrac{2}{n}\bar x\sum x_i + \bar x^2\\[6pt] &= \dfrac1n\sum x_i^2 - 2\bar x^2 + \bar x^2\\[6pt] &= \dfrac1n\sum x_i^2 - \bar x^2. \end{array} \]

So with \(S_1=\sum x_i\) and \(S_2=\sum x_i^2\),

\[ \sigma^2 = \frac{S_2}{n} - \left(\frac{S_1}{n}\right)^2. \]

\(S_1\), \(S_2\) and \(n\) are ordinary sums, so a group's summary can be kept in exactly three columns:

CREATE TABLE t2(
  group_value INT UNSIGNED NOT NULL PRIMARY KEY,
  group_count INT UNSIGNED NOT NULL DEFAULT 0,
  group_sum DOUBLE NOT NULL DEFAULT 0,
  group_sum2 DOUBLE NOT NULL DEFAULT 0
);

Feeding it one row at a time is a plain accumulation, with no reference to any other row in the group:

INSERT INTO t2 (group_value, group_count, group_sum, group_sum2)
SELECT group_value, 1, value, value * value
FROM t1
AS new
ON DUPLICATE KEY UPDATE
  group_count = t2.group_count + 1,
  group_sum   = t2.group_sum + new.group_sum,
  group_sum2  = t2.group_sum2 + new.group_sum2;

SELECT group_value GVal, group_count GCnt, group_sum / group_count GAvg,
       SQRT(group_sum2 / group_count - POW(group_sum / group_count, 2)) GStd
FROM t2;

Because \(S_1\) and \(S_2\) are additive, regrouping the aggregated table further needs nothing beyond another SUM(). Splitting the five groups into (group_value <= 2) and its complement, for instance:

SELECT group_value <= 2 GVal, SUM(group_count) GCnt,
       SUM(group_sum) / SUM(group_count) GAvg,
       SQRT(SUM(group_sum2) / SUM(group_count) - POW(SUM(group_sum) / SUM(group_count), 2)) GStd
FROM t2
GROUP BY 1;

reproduces exactly what STD() would compute over the matching rows of \(t_1\) - no join, no special-case formula, because every quantity involved is still a plain sum of sums.

Welford's Online Algorithm

The sum-of-squares approach is simple, but it has a numerical problem: \(S_2/n\) and \(\bar x^2\) can both be huge, nearly equal DOUBLE values whenever the mean is large compared to the spread of the data, and subtracting two close, large numbers cancels most of their significant digits. For values around \(10^8\) with a spread of only a few units, \(S_2/n\) and \(\bar x^2\) agree in about 14 digits before the subtraction, leaving almost nothing of the true variance - rounding error alone can push the result slightly negative, which turns SQRT() into NULL.

Welford's online algorithm [Welford] avoids squaring the raw values altogether by tracking the running mean and the running sum of squared deviations from that mean instead of the two sums \(S_1\), \(S_2\). The chapter on rolling variance derives its update rule in full from the definition of variance; only its final normalization changes here, since STD() divides by \(n\), not \(n-1\). Writing \(\bar x_n\) for the running mean and \(M_{2,n}:=\sum_{i=1}^n (x_i-\bar x_n)^2\) for the running sum of squared deviations, the update after a new value \(x_n\) arrives is

\[ \bar x_n = \bar x_{n-1} + \frac{x_n-\bar x_{n-1}}{n}, \qquad M_{2,n} = M_{2,n-1} + (x_n-\bar x_{n-1})(x_n-\bar x_n), \]

and the population variance is \(\sigma_n^2=M_{2,n}/n\) exactly - only the divisor differs from the sample form \(s_n^2=M_{2,n}/(n-1)\) derived in that chapter. As pseudocode, using only the value that just arrived, the mean before the update, and the mean after it:

n = n + 1
delta = x - mean
mean = mean + delta / n
delta2 = x - mean
m2 = m2 + delta * delta2

variance = m2 / n
std = sqrt(variance)

A group's summary now needs the running mean and \(M_2\) instead of the two sums:

CREATE TABLE t3(
  group_value INT UNSIGNED NOT NULL PRIMARY KEY,
  group_count INT UNSIGNED NOT NULL DEFAULT 0,
  group_mean DOUBLE NOT NULL DEFAULT 0,
  group_m2 DOUBLE NOT NULL DEFAULT 0
);

Translating the update rule to SQL needs one extra step: group_mean is only overwritten once, so the expression for group_m2 has to compute the new mean inline instead of referring to the (not yet updated) group_mean column:

INSERT INTO t3 (group_value, group_count, group_mean, group_m2)
SELECT group_value, 1, value, 0
FROM t1
AS new
ON DUPLICATE KEY UPDATE
  group_count = t3.group_count + 1,
  group_m2    = t3.group_m2 + (new.group_mean - t3.group_mean)
                * (new.group_mean - (t3.group_mean + (new.group_mean - t3.group_mean) / group_count)),
  group_mean  = t3.group_mean + (new.group_mean - t3.group_mean) / group_count;

SELECT group_value GVal, group_count GCnt, group_mean GAvg, SQRT(group_m2 / group_count) GStd
FROM t3;

This reproduces the same AVG() / STD() baseline as the sum-of-squares version, just without ever squaring a raw value.

Combining Two Aggregated Groups

Regrouping \(t_3\) is not a plain SUM() the way it was for \(t_2\), because \(M_2\) is only the sum of squared deviations from its own group's mean - two groups with different means cannot simply have their \(M_2\) added together. Combining group \(X\) (\(n_x\) values, mean \(\bar x\), sum of squared deviations \(M_{2,x}\)) with group \(Y\) (\(n_y\), \(\bar y\), \(M_{2,y}\)) into \(Z=X\cup Y\) needs its own formula.

The combined count and mean follow directly from summing every value in \(Z\):

\[ n_z = n_x+n_y,\qquad \bar z = \frac{n_x\bar x+n_y\bar y}{n_z}. \]

For \(M_{2,z}=\sum_{X}(x_i-\bar z)^2+\sum_Y(y_i-\bar z)^2\), the same trick as above - inserting and subtracting each group's own mean - applies to both sums:

\[ \sum_X(x_i-\bar z)^2 = \sum_X\big[(x_i-\bar x)+(\bar x-\bar z)\big]^2 = M_{2,x} + n_x(\bar x-\bar z)^2, \]

since \(\sum_X(x_i-\bar x)=0\) removes the cross term again, and likewise \(\sum_Y(y_i-\bar z)^2=M_{2,y}+n_y(\bar y-\bar z)^2\). Adding both halves gives the pairwise combination formula [Chan]:

\[ M_{2,z} = M_{2,x} + n_x(\bar x-\bar z)^2 + M_{2,y} + n_y(\bar y-\bar z)^2, \qquad \sigma_z^2 = \frac{M_{2,z}}{n_z}. \]

In SQL, \(\bar z\) has to be computed first, since every other term depends on it:

SELECT g.GVal,
       SQRT(SUM(t3.group_m2 + t3.group_count * POW(t3.group_mean - g.GMean, 2)) / SUM(t3.group_count)) GStd
FROM t3
JOIN (
  SELECT group_value <= 2 GVal,
         SUM(group_mean * group_count) / SUM(group_count) GMean
  FROM t3
  GROUP BY 1
) g ON (t3.group_value <= 2) = g.GVal
GROUP BY g.GVal;

which reproduces the same regrouped STD() baseline as the plain-sum version above, at the cost of the extra join that \(t_2\)'s additive statistics never needed.

Common Pitfalls

• STD(), STDDEV() and VARIANCE() are population statistics (divide by \(n\)); use STDDEV_SAMP() / VAR_SAMP() if the sample form (divide by \(n-1\)) is what's actually needed, and adjust the final division in every query above accordingly.
• The sum-of-squares method is only safe when the values are small enough, or spread out enough, that \(S_2/n\) and \(\bar x^2\) don't cancel each other down to rounding noise; Welford's running-mean form should be preferred whenever that isn't guaranteed.
• VALUES(col_name) inside ON DUPLICATE KEY UPDATE is deprecated since MySQL 8.0.20; aliasing the inserted row with AS new and referring to new.col_name, as done throughout this article, is the current replacement.
• Regrouping \(t_2\) or \(t_3\) further only works along a key that is a deterministic function of columns already stored in the aggregate table, such as group_value <= 2 here; regrouping by anything not derivable from that key requires re-aggregating from the raw data.
• Both tables assume rows are only ever inserted, never edited or deleted; removing or changing an individual value afterwards leaves the running aggregate silently wrong.

References

• Welford B. P. Welford (1962) Note on a Method for Calculating Corrected Sums of Squares and Products, Technometrics, 4(3)
• Chan T. F. Chan, G. H. Golub, R. J. LeVeque (1979) Updating Formulae and a Pairwise Algorithm for Computing Sample Variances, Technical Report STAN-CS-79-773, Stanford University