root/trunk/mathextra/leastsqr.d

// Cholesky decomposition
module mathextra.leastsqr;
import std.math;

debug private import std.stdio;

unittest {
    debug writefln("--- UnitTest: " __FILE__ " ---");
} 


/**
 *  Perform a multiple linear regression (least squares fit).
 *
 *
  Given a linear equation with unknown coefficients a[],
     a[0]*f[0](i) + a[1]*f[1](i) + ... = yfunc(i)
  this function uses a Cholesky decomposition to find the values 
  of a[] which minimize the
  sum of squares, defined as:
  
  sumsq = $(BIGSUM i) $(POWER ( y[i] - $(BIGSUM j) a[j] * f[j](i) ), 2)
  
  The f[]() and yfunc() functions must be linearly independent.
  
  
  If the Cholesky decomposition is not positive definite,
  the corresponding a[] entry will be nan; this indicates that
  there was insufficient information to determine that value.
 *
 * Params:
 *  N       Number of coefficients.
 *  numpoints Number of points in yfunc and xfunc.
 *  f       Array of delegates. f[j](i) is the i'th data point for coefficient j.
 *  yfunc   yfunc(i) is the i'th data point 
 *
 * Returns:
 *          The sum of squares, sumsq. The coefficients are returned
 *          in solution.
 */
template leastSquaresFit(int N)
{
    real leastSquaresFit(real solution[N], int numpoints, real delegate (int) yfunc, real delegate (int) [N] f ...) {
        return Regression!(N).fitCholesky(solution, numpoints, yfunc, f);
    }
}

template Regression(int N) {

/*
  Given a positive-definite matrix A,
 constructs its cholesky decomposition,
 A = L. L<sup>T</sup>
 
 Only the upper triangular part of A is required.
 eg A[2][3].
 Returns result in lower triangle of A. Diagonals
 are returned in diag.
*/
// diag[0] will be 0 for any row which is not positive definite.
// The corresponding row & column in the matrix will also be all zeros.
void CholeskyDecompose(real a[N][N], real diag[N])
{
    for (int i=0; i<N; ++i) {
        real sum = a[i][i];
        for (int k=i-1; k>=0; --k) {
            sum -= a[i][k] * a[i][k];
        }
        if (sum >0.0) {
            diag[i]=sqrt(sum);
           
            for (int j=i+1; j<N; ++j) {
                sum=a[i][j];
                for (int k=i-1; k>=0; --k) {
                    sum -= a[i][k] * a[j][k];
                }
                a[j][i]=sum/diag[i];
            }
        } else {
            // not positive definite (could be caused by rounding errors)
            diag[i] = 0;
            // make this whole row zero so they have no further effect
            for (int j=i+1; j<N; ++j) a[j][i] = 0;
        }
    }
}

void CholeskySolve(real a[N][N], real diag[N], real b[N], real x[N])
{
    for (int i=0; i<N; ++i) {
        if (diag[i]>0)  {
            real sum=b[i];
            for (int k=i-1; k>=0; k--) sum-=a[i][k]*x[k];
            x[i] = sum/diag[i];
        } else x[i] = 0; // skip pos definite rows
    }
    for (int i=N-1; i>=0; --i) {
        if (diag[i]>0) {
            real sum=x[i];
            for (int k=i+1; k<N; ++k) sum-=a[k][i]*x[k];
            x[i] = sum/diag[i];
        } else x[i]=0; // skip pos definite rows
    }
    // Convert failed columns in solution to NANs if required.
    for (int i=0; i<N; ++i) {
        if (diag[i]!>0) x[i]=real.nan;
    }
}

// The x and y values are calculated from indices.
// If for any n, the functions return nan, that point will be ignored.
real fitCholesky(real solution[N], int numpoints, real delegate (int) yfunc, real delegate (int) [N] f ...)
{
    real [N][N] m;
    real [N] b;
    // Clear the matrices
    for (int i=0; i<N; ++i) {
        m[i][i..$] = 0;
    }
    b[0..$] = 0.0;
    
    real sumsq_y = 0;

    for (int q=0; q<numpoints; ++q) {
        real [N] v;
        real y = yfunc(q);
        real anynans = y;
        for (int i=0; i<N; ++i) {           
            v[i] = f[i](q);
            anynans *= v[i];
        }
        if (isnan(anynans)) continue;
        sumsq_y += y*y;
        for (int i=0; i<N; ++i) {
            b[i] += y * v[i];
            for (int j=i; j<N; ++j) {
                 m[i][j] += v[i]*v[j];
            }
        }
    }
    real diag[N];
    CholeskyDecompose(m, diag);
    CholeskySolve(m, diag, b, solution);
    
    real sumsq=0;
    for (int i=0; i<N; ++i) {
        if (isnan(solution[i])) continue;
        for (int j=i+1; j<N; ++j) {
            if (isnan(solution[j])) continue;
            sumsq += m[i][j]  * solution[i]*solution[j];
        }
    }
    sumsq*=2; // account for the lower triangular part
    for (int i=0; i<N; ++i) {
        if (isnan(solution[i])) continue;
        sumsq += m[i][i]*solution[i]*solution[i];
        sumsq -= 2 * b[i] * solution[i];
    }
    return sumsq + sumsq_y;
 }
 
/+
// returns the sum of squares. Used only for unit tests.
real evaluateFit(real solution[N], int numpoints, real delegate (int) yfunc, real delegate (int) [N] xfunc ...)
{
    static real sqr(real x) { return x*x; }
       
    real sumsq=0;
    for (int i=0; i<numpoints; ++i) {
        real a=0;
        for (int j=0; j<N; ++j) {
            if (!isnan(solution[j]))  a += solution[j] * xfunc[j](i);
        }
        sumsq += sqr( a - yfunc(i) );
    }
    return sumsq;
}
+/


} // end Regression


unittest {
    real x[];
    real y[];
    x.length = 500;
    y.length = 500;
    real intercept = 580; 
   for (int i=0; i<x.length; ++i) {
     x[i] = i*12.82344 + 458743.3;
     y[i] = x[i]*x[i]*PI + intercept;
   }

   // fit y = a[0]x^2 + a[1]*1
   
   real getx(int n) { return x[n]*x[n]; }
   real gety(int n) { return y[n]; }
   
   real [2] a;
   
   real delegate (int)[] delegs = void;
   leastSquaresFit!(2)(a, x.length, &gety, &getx, delegate real (int n) { return 1; });
 
  assert(feqrel(a[0],PI)>=double.mant_dig-5);
  assert(fabs(a[1]-intercept) < 0.002);
}