fast & efficient least squares fit algorithm in C?

Question

I am trying to implement a linear least squares fit onto 2 arrays of data: time vs amplitude. The only technique I know so far is to test all of the possible m and b points in (y = m*x+b) and then find out which combination fits my data best so that it has the least error. However, I think iterating so many combinations is sometimes useless because it tests out everything. Are there any techniques to speed up the process that I don't know about? Thanks.

Answer 1

Try this code. It fits y = mx + b to your (x,y) data.

The arguments to linreg are

linreg(int n, REAL x[], REAL y[], REAL* b, REAL* m, REAL* r)  n = number of data points x,y  = arrays of data *b = output intercept *m  = output slope *r = output correlation coefficient (can be NULL if you don't want it) 

The return value is 0 on success, !=0 on failure.

Here's the code

#include "linreg.h" #include <stdlib.h> #include <math.h>                           /* math functions */  //#define REAL float #define REAL double   inline static REAL sqr(REAL x) {     return x*x; }   int linreg(int n, const REAL x[], const REAL y[], REAL* m, REAL* b, REAL* r){     REAL   sumx = 0.0;                      /* sum of x     */     REAL   sumx2 = 0.0;                     /* sum of x**2  */     REAL   sumxy = 0.0;                     /* sum of x * y */     REAL   sumy = 0.0;                      /* sum of y     */     REAL   sumy2 = 0.0;                     /* sum of y**2  */      for (int i=0;i<n;i++){          sumx  += x[i];                sumx2 += sqr(x[i]);           sumxy += x[i] * y[i];         sumy  += y[i];               sumy2 += sqr(y[i]);      }       REAL denom = (n * sumx2 - sqr(sumx));     if (denom == 0) {         // singular matrix. can't solve the problem.         *m = 0;         *b = 0;         if (r) *r = 0;             return 1;     }      *m = (n * sumxy  -  sumx * sumy) / denom;     *b = (sumy * sumx2  -  sumx * sumxy) / denom;     if (r!=NULL) {         *r = (sumxy - sumx * sumy / n) /    /* compute correlation coeff */               sqrt((sumx2 - sqr(sumx)/n) *               (sumy2 - sqr(sumy)/n));     }      return 0;  } 

Example

You can run this example online.

int main() {     int n = 6;     REAL x[6]= {1, 2, 4,  5,  10, 20};     REAL y[6]= {4, 6, 12, 15, 34, 68};      REAL m,b,r;     linreg(n,x,y,&m,&b,&r);     printf("m=%g b=%g r=%g\n",m,b,r);     return 0; } 

Here is the output

m=3.43651 b=-0.888889 r=0.999192     

Here is the Excel plot and linear fit (for verification).

All values agree exactly with the C code above (note C code returns r while Excel returns R**2).

Answer 2

There are efficient algorithms for least-squares fitting; see Wikipedia for details. There are also libraries that implement the algorithms for you, likely more efficiently than a naive implementation would do; the GNU Scientific Library is one example, but there are others under more lenient licenses as well.