int lsq_zero (n, mat, vec) int lsq_accum (n, mat, vec, a, w) int lsq_solve (n, mat, vec, sol) int lsq_cfill (n, mat, c, sol) int n, c; real mat[n*n], vec[n], sol[n], a[n+1], w; DescriptionThese routines provide a low level interface to solving linear least squares problems using the Normal Equations (see e.g. Numerical Recepies, Press et al.). These routines have the advantage that they don’t need the extra temporary arrays containing the data, i.e. the Design Matrix, and directly accumulate the data in the appropriate spots in the symmetric Normal Matrix. lsq_accum accumulates the data into the normal equations matrix that is inverted in lsq_solve. Before the data is accumulated, the matrix and right hand side vector need to be reset using lsq_zero. During accumulation, a weight, w, can be given to each datapoint (Note: this is a linear weight, by which the data is multiplied. If you’re used to enter the dispersion, you should enter 1/sigma^2). On output the matrix mat contains the inverse of the normal matrix, and hence its diagonal elements the square of the errors of the fitted parameters. The fitted parameters themselves are returned in the array sol. ExampleIn this example a large 2D image matrix is fitted with an intensity gradient of the form I(x,y)=a+bx+cy: double mat[3*3],vec[3],a[4], real data[nx*ny]; lsq_zero(3,mat,vec); for (ix=0; ix<nx; ix++) for (iy=0; iy<ny; iy++) { a[0] = 1.0 a[1] = ix; a[2] = iy; a[3] = data[ix*ny+iy]; lsq_accum(3,mat,vec,a,1.0); } lsq_solve(3,mat,vec,a); printf("a+bx+cy: a=%f b=%f c=%f\n", a[0],a[1],a[2]);

~/src/kernel/misclsq.c

29-sep-90createdPJT 19-feb-92updated doc, and properly redfined the weightsPJT