numerically stable covariance algorithm

covariance.c

int covariance(int n, int m, double data[], int strides[2], char mode, double matrix[])
/* This algorithm is described in:
 * B.P. Welford:
 * "Note on a method for calculating corrected sums of squares and products."
 * Technometrics 4(3): 419-420 (1962).
 * Also see:
 * Peter M. Neely:
 * "Comparison of several algorithms for computation of means, standard
 * deviations and correlation coefficients."
 * Communications of the ACM 9(7): 496-499 (1966).
 */
{ int i, j, k;
  double* p;
  double* q;
  double* q1;
  double* q2;
  double x1;
  double x2;
  const int stride1 = strides[0];
  const int stride2 = strides[1];
  const int denominator = (mode=='u') ? n-1 : n;
  /* 'u': Unbiased estimate; otherwise maximum-likelihood estimate */
  double* average = malloc(m*sizeof(double));
  if(!average) return 0;

  /* Initialize to zero */
  p = matrix;
  for (i = 0; i < m; i++)
  { average[i] = 0.0;
    for (j = 0; j < m; j++, p++) *p = 0.0;
  }

  /* Calculate the sums of squares */
  for (k = 0; k < n; k++)
  { const double scale = k+1.0;
    const double factor = k/scale;
    q = data + k*stride1;
    q1 = q;
    for (i = 0; i < m; i++, q1+=stride2)
    { p = matrix + i*m;
      x1 = (*q1) - average[i];
      q2 = q;
      for (j = 0; j <= i; j++, p++, q2+=stride2)
      { x2 = (*q2) - average[j];
        *p += factor*x1*x2;
      }
    }
    for (i = 0; i < m; i++, q+=stride2)
      average[i] = factor * average[i] + (*q) / scale;
  }
  free(average);
 
  /* Scale, and copy to the upper half of the matrix */
  for (i = 0; i < m; i++)
  { p = matrix + i*m;
    q = matrix + i + (i+1)*m;
    p[i] /= denominator;
    for (j = i+1; j < m; j++, q+=m)
    { (*q) /= denominator;
      p[j] = (*q);
    }
  }
  return 1;
}