11
\$\begingroup\$

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;
}
\$\endgroup\$
2
  • 2
    \$\begingroup\$ Have you profiled or do a benchmark to find out where in your code most of the time is spent? In most programs involving I/O, the I/O is the bottleneck. If you can use block reading. The getline doesn't block read; it actually is scanning for a ending sentinel; which slows down the input process. \$\endgroup\$ Commented Jul 10, 2023 at 23:55
  • 3
    \$\begingroup\$ I have rolled back your last edit. Please do not update the code in your question to incorporate feedback from answers, doing so goes against the Question + Answer style of Code Review. This is not a forum where you should keep the most updated version in your question. Please see what you may and may not do after receiving answers. \$\endgroup\$
    – Heslacher
    Commented Jul 11, 2023 at 8:54

5 Answers 5

10
\$\begingroup\$

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.

\$\endgroup\$
6
\$\begingroup\$

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.

\$\endgroup\$
2
  • \$\begingroup\$ Regarding item count: on some systems (e. g. Windows), size_t and/or unsigned long cannot hold the amount of records that may occur in a file. Possible alternatives are off_t (the return type of ftello() which is meant to supersede ftell()) and uintmax_t which are guaranteed to hold the maximum supported file size on a platform. \$\endgroup\$ Commented Jul 11, 2023 at 15:19
  • \$\begingroup\$ @DavidFoerster off_t, ftello() are not part of the standard C library. Concerning "on some systems (e. g. Windows)", I'd say the issue is a compiler one, not an OS one. I find no specified guarantee that uintmax_t is sufficiently large for all file sizes, yet it is very reasonable and practical to assume that. \$\endgroup\$ Commented Jul 11, 2023 at 15:32
6
\$\begingroup\$

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.

\$\endgroup\$
3
  • 1
    \$\begingroup\$ Hmmm, with no return on main() then "function returns a value of 0." § 5.1.2.2.3 1. EXIT_SUCCESS is not defined as 0 but with use with exit(): "EXIT_SUCCESS which expand to integer constant expressions that can be used as the argument to the exit function" and implies it might differ from 0: "If the value of status is zero or EXIT_SUCCESS": Curiously, the spec does not talk about return EXIT_SUCCESS; in main(). I wonder 1) if that is UB 2) Has any implementation done EXIT_SUCCESS != 0? Think I will stay with return 0; for now. \$\endgroup\$ Commented Jul 12, 2023 at 6:44
  • 1
    \$\begingroup\$ Ooh, that's worthy of a language-lawyer question over on Stack Overflow, I think. \$\endgroup\$ Commented Jul 18, 2023 at 7:12
  • 1
    \$\begingroup\$ The spec does define the behaviour of return EXIT_SUCCESS; albeit indirectly: "a return from the initial call to the main function is equivalent to calling the exit function with the value returned by the main function as its argument" (my emphasis). I guess that EXIT_SUCCESS is to successful exit somewhat as NULL is to null pointers: a macro that doesn't necessarily expand to 0 but causes the same behaviour as 0. \$\endgroup\$ Commented Jul 18, 2023 at 7:29
3
\$\begingroup\$

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.

\$\endgroup\$
1
  • \$\begingroup\$ Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on Code Review Meta, or in Code Review Chat. Comments continuing discussion may be removed. \$\endgroup\$
    – Malachi
    Commented Jul 17, 2023 at 16:39
0
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.