Is Highams' computation of mean worth the price?

Question

In Accuracy and Stability of Numerical Algorithms, equation 1.6a, Higham gives the following update formula for the mean: $$ M_{1} := x_1, \quad M_{k+1} := M_{k} + \frac{x_k - M_k}{k} $$ Ok, one benefit of this update is that it doesn't overflow, which is a potential problem when we naively compute $\sum_{i=1}^{k} x_k$ and then divide by $k$ before any use.

But the naive sum requires $N-1$ adds and 1 division, and potentially vectorizes very nicely, whereas the update proposed by Higham requires $2N-2$ additions and $N-1$ divisions, and I'm unaware of any assembly instruction that can vectorize all this.

So it Higham's update formula worth using? Are there benefits I'm not seeing?

Note: Higham gives the following generic advice (Section 1.18 "Designing Stable Algorithms):

"It is advantageous to express update formulae as newvalue = oldvalue + smallcorrection if the small correction can be computed with many correct significant figures."

Update 1.6a does take this form, but it's not clear to me that the small correction can be given to many significant figures.

Edit: I found an empirical study of the performance of various methods of computing means, and 1.6a came highly recommended; see

Robert F. Ling (1974) Comparison of Several Algorithms for Computing Sample Means and Variances, Journal of the American Statistical Association, 69:348, 859-866, DOI: 10.1080/01621459.1974.10480219

However, it still wasn't clear to me after reading that paper that the update was worth the price; and in any case, I was hoping for a worst and average case bound by accumulating rounding errors.

Answer 1

Higham's algorithm seems very useful for online algorithms or algorithms in which you have limited storage capacity such as edge processing. In these situations it is probably always worth the cost.

But, to address somewhat your question as posed, I implemented an SSE2 version which I think captures your question:

#include <chrono>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <xmmintrin.h>

int main(){
  const int num_count = 100000;

  alignas(16) float data[num_count];
  for(int i=0;i<num_count;i++)
    data[i] = rand()%100000+rand()/(double)RAND_MAX;


  {
    const auto t1 = std::chrono::high_resolution_clock::now();
    float sum=0;
    for(int i=0;i<num_count;i++){
      sum += data[i];
    }
    const float mean=sum/num_count;

    const auto t2 = std::chrono::high_resolution_clock::now();
    const auto time_span1 = std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1);
    std::cout << "Exec time:\t" << time_span1.count() << " s\n";
    std::cout << "Mean = "<<std::setprecision(20)<<mean<<std::endl;
  }

  {
    const auto t1 = std::chrono::high_resolution_clock::now();
    __m128 mean = _mm_load_ps(&data[0]);
    for (int i=4;i<num_count;i+=4){
      const __m128 x    = _mm_load_ps(&data[i]);
      const __m128 diff = _mm_sub_ps(x,mean);
      const __m128 k    = _mm_set1_ps(i/4);
      const __m128 div  = _mm_div_ps(diff, k);
      mean              = _mm_add_ps(mean, div);
    }
    float result[4];
    _mm_store_ps(result, mean);
    const float tmean = (result[0] + result[1] + result[2] + result[3])/4; //I'm suspicious about this step: probably throws away precision
    const auto t2 = std::chrono::high_resolution_clock::now();
    const auto time_span1 = std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1);
    std::cout << "Exec time:\t" << time_span1.count() << " s\n";
    std::cout << "Mean = "<<std::setprecision(20)<<tmean<<std::endl;
  }


}

}

and observe that

Exec time:  0.000225851 s
Simple Mean = 49891.23046875
Exec time:  0.0003759360000000000002 s
Higham Mean = 49890.26171875

Higham's mean takes a little longer to calculate and the values differ by what might be a non-negligible amount, though the accuracy you need really depends on your implementation.

Answer 2

Suppose you have $m$ numbers and the vector length $l$ divides $n$, i.e., $m=ln$. The you can apply Higham's scheme to each group of $n$ numbers. This vectorizes nicely, especially if the numbers are interleaved, i.e., all the first elements are contiguous in memory, followed by all the second elements, etc. When you have computed the average of each group you have $l$ numbers left. You apply a sequential implementation of Higham's scheme and then you are done. If the vector length does not divide $n$, then pad the data with zeros, split the data into $l$ groups with $n$ numbers each and proceed as before. When you have the average $\mu_j$ of each group, you compute $(n/m)\mu_j$ and simply add these numbers using the naive summation algorithm.

---

Whether Higham's formula is worth using depends on the situation. When overflow is an issue, then the naive formula is useless. Both formulas have low arithmetic intensity, and I expect Higham's formula to run very slowly simply because the division has very high latency. Personally, I find that speed is very secondary compared with reliability and accuracy.

---

Higham's primary goal in this very first chapter is simply to alert the reader to the fact that computing with floating point numbers is profoundly different from real arithmetic. His advice about corrections is quoted correctly, but I would have written it a bit differently. Specifically, it does not matter if the correction is computed with a large relative error as long as the correction is small relative to the original value.