Calculating standard deviation

Question

I am getting about 215 MB/s throughput on this implementation of Welford's algorithm. My hard drive is rated for 400 MB/s reads, so I am wondering how I can optimise this. I assume the bottleneck is in getline or atof as 400 MB/s throughput gives my 4 GHz processor 640 cycles per double to work with, which should be more than enough for the computational part of the loop. Are there faster alternatives?

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc, char **argv)
{
    if (argc != 1) {
        printf("Usage: %s < data\n", argv[0]);
        return EXIT_FAILURE;
    }

    size_t count = 0;
    double mean = 0;
    double m2 = 0;

    char *line = NULL;
    size_t len = 0;
    ssize_t nread;

    while ((nread = getline(&line, &len, stdin)) != -1) {
        double newValue = atof(line);
        count++;
        double delta = newValue - mean;
        mean += delta / count;
        m2 += delta * (newValue - mean);
    }

    printf("%lf\n%lf\n", mean, sqrt(m2 / count));

    free(line);

    return 0;
}

Answer 1

Avoid intermediate copies

By using getline() you make a copy of some characters into a temporary buffer, then you call atof() on that. There is some overhead in that. You could consider doing this instead:

while(scanf("%lf ", &newValue) == 1) {
    …
}

Which may or may not be faster. Note that getline() also allocates heap memory for the buffer, which probably doesn't matter much since it will only do it a few times at most, but it could be avoided by using a fixed-size array; after all, numbers are only so big.

Use a better floating point number parser

The standard library might not have an optimal algorithm to convert strings to floating point values, and apart from that it might take into account locales, which adds some overhead. In recent years some advances have been made in string-to-float algorithms, and some (standard) libraries also have functions to parse numbers without bothering with locales. The C++ library for example has gotten std::from_chars(). You might also want to look at external libaries, like Daniel Lemire's fast_double_parser.

Use multiple threads

Even if you can make the parsing so fast that your program is purely I/O bound, you might still get a benefit from making it multi-threaded. That's because SSDs nowadays can work faster if they receive multiple I/O requests to different parts of the SSD.

Missing error checking

Your should check after the while-loop that feof(stdin) == true, to verify that you actually read all the way to the end of the file. If it's false, that means an error occured, in which case you should not print out an incorrect standard deviation, but print an error message to stderr and return EXIT_FAILURE.

Note that parsing a number might fail as well (maybe the file containing the data was corrupted for example, or it's not a file with numbers at all), so it would be better if you could check for those kinds of errors as well, although it will come at some cost to performance.

Answer 2

Missing assessment

OP's goal is "how I can optimise this.", yet post lacks a test harness. Posting code that reports the performance allows us all to post real results rather than a bunch of theoretical improvements and ideas.

Fixed point printing

printf("%lf\n%lf\n", mean, sqrt(m2 / count)); hints that precision is not the real the goal.

Printing floating point numbers with fixed point format prints about 25% of double as 0.000000 and 25% of large values with excessive precision.

Use "%.16e", "%.17g" or "%a" for more informative results.

---

Simplify algorithm

Consider base-lining with the basic naive algorithm.

while ((nread = getline(&line, &len, stdin)) != -1) {
    double newValue = atof(line);
    sum += newValue;
    sum_squared += newValue * newValue;
    count++;
}
printf("Average: %g\n", sum/count);
printf("Variance: %g\n", (sum_squared − (sum*sum) / count) / (count − 1));

Depending on the data set, the 64-bit double math with the simpler algorithm may render a sufficiently precise and speedier answer.

Part of the reason for the 15+ decimal digits of precision for double is so one does not need to perform precision preservating code and instead go for direct fast encoding.

float

Using float vs. double is often faster (4x even). Much depends on the underlying hardware, yet worth a try.

Time for I/O

Timing a simple read-only would help gauge how much time is spent reading.

clock_t c0 = clock();
while ((nread = getline(&line, &len, stdin)) != -1) {
  ;
}
clock_t c1 = clock();
printf("time: %g seconds\n", 1.0*(c1 - c0)/CLOCKS_PER_SEC);

---

size_t?

Unclear why code uses size_t to count the values. size_t correlates with memory size. I suppose some type that goes along with file size is useful. unsigned long is a reasonable alternative given ftell() returns long.

With the naïve approach, might as well use uintmax_t.

Answer 3

The calculations appear to use the numerically stable formulae from equations (4) and (24) in Incremental calculation of weighted mean and variance - well done for avoiding the mathematician's version that suffers from catastrophic cancellation in floating-point arithmetic.

However, the correspondence isn't immediately obvious here:

        m2 += delta * (newValue - mean);

Because we assigned delta == newValue - mean, it looks like this is the same as delta * delta. However, because mean has been modified in between, it's not equivalent. Perhaps it's worth renaming delta to delta_pre and creating an new delta_post to make this absolutely clear (and prevent a subsequent maintainer from "correcting" the code):

        double new_value = atof(line);
        double const delta_pre = new_value - mean;
        mean += delta_pre / ++count;

        double const delta_post = new_value - mean;
        m2 += delta_pre * delta_post;

---

Minor things:

• Error message should go to stderr, not stdout:

      if (argc != 1) {
          fprintf(stderr, "Usage: %s < data\n", argv[0]);
  

  The message is perhaps misleading as there are other ways of providing data on standard input, such as here-strings or pipes.

• We could return EXIT_SUCCESS instead of magic number 0 at the end, to be more consistent with the error return. Or just omit the return statement - main() is magic and will return success just by executing off the end.

Answer 4

Profile to Find the Bottleneck

Others have already made this point, but you need to find where your program is spending all its time, and focus on that.

Look for Ways to Speed up the I/O

This is going to be implementation-dependent, but in many cases, fscanf(), which doesn’t need to copy or heap-allocate, can be faster than getline(). In fact, you currently allocate and free a buffer from the heap for every data point.

Not in the standard library, your OS, or a third-party library like glib, probably lets you memory-map a regular file. That lets you scan for tokens and pass the beginning and end of each token to strtod. This can be extremely fast because there is no buffering or copying, just mapping the pages of the file on disk to memory similarly to the virtual memory manager. (G. Sliepen had an excellent series of comments, now removed, pointing out the ways that might not be true.)

Store the Data Points in a Dynamic Array

Instead of reading each value in one at a time, you want to store all the scanned values in memory (maximizing the I/O throughput) and then process them at once (allowing parallelization and SIMD).

It’s possible to implement a resizing vector of double in C, and that’s probably the data structure you want for this, although you might also make two passes through the file in memory, one to count the items, and then another to fill the array.

Parallelizing the Calculation

The version of the calculation that you’re using is incremental, and therefore inherently serial. If you have all the data in memory, however, there is no need to use an incremental formula. You can make two passes, one reduction to calculate the mean, and then another to calculate the variance as a map-reduction (the sum of (xᵢ-μ)² for all xᵢ), and finally calculate σ from the variance. You will then be adding the squares of small differences, which all have the same sign, and not get any subtractive cancellation from this step (although there could be large data points with opposite signs making the calculation of the mean numerically unstable). This version of the calculation sacrifices some stability for being able to use SIMD.

You can additionally split the both passes between threads, with OpenMP, or calculate the variance of each subrange of the input and then the total variance from those. Since the mean and variance can be calculated from simple sums of terms, you could use the reduction operation of OpenMP for this (or in the C++ standard library, std::transform_reduce with a parallel execution policy), rather than a weighted-mean and weighted-variance formula. Otherwise, you would want to partition the input array into one slice per worker thread.

While the way I outlined above sacrifices some numerical stability for the speed of SIMD, you might figure out an incremental formula that adds four terms at once, or however many your hardware can efficiently vectorize.

Answer 5

The standard library functions often have features which make them nice to use but also introduce overhead. If this is for competitive programming, your input number format is probably very simple, and you can take advantage of this by handling everything manually. Allocate a buffer, say a megabyte, do an fread into it, and parse the buffer with simple logic like if (buffer[position] >= '0' && buffer[position] <= '9'). Make sure to handle the case of a number spanning the end of one read and beginning of another.

If this is going to run in production, I'd keep it as it is, except I'd replace atof with strtod in order to do a bit of error checking. It's slower than manually parsing a buffer, but it's harder to make a programming mistake this way.