Table of Contents

Name

lsq_zero, lsq_accum, lsq_solve - least squares fitting utilities

Synopsis


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]);

See Also

lsqfit(3NEMO) , matinv(3NEMO) , mpfit(3NEMO)

Author

Peter Teuben

Files


~/src/kernel/misc      lsq.c

Update History


29-sep-90    created      PJT
19-feb-92    updated doc, and properly redfined the weights    PJT


Table of Contents