rebase fitter on balance and Jacobian instead of gradient and Hessian

Author: JohnMount

Count/src/com/mzlabs/count/ctab/CTab.java

@@ -15,7 +15,7 @@
 import com.mzlabs.count.op.impl.ThreadedSum;
 import com.mzlabs.count.op.iter.OrderStepperTot;
 import com.mzlabs.count.zeroone.ZeroOneCounter;
-import com.mzlabs.fit.GLMModel;
+import com.mzlabs.fit.DirectPoissonJacobian;
 import com.mzlabs.fit.NewtonFitter;
 
 
@@ -216,7 +216,7 @@ public static void main(final String[] args) {
 		for(int n=1;n<=9;++n) {
 			final CTab ctab = new CTab(n,true);
 			//final NewtonFitter lf = new NewtonFitter(new SquareLossOfExp());
-			final NewtonFitter lf = new NewtonFitter(GLMModel.PoissonLink);
+			final NewtonFitter lf = new NewtonFitter(DirectPoissonJacobian.poissonLink);
 			final int tLast = (n*n-3*n+2)/2;
 			for(int total=0;total<=tLast;++total) {
 				final Date startTime = new Date();

Count/src/com/mzlabs/fit/BalanceBasedJacobian.java

@@ -0,0 +1,36 @@
+package com.mzlabs.fit;
+
+
+/**
+ * This returns a set of equations that should be zero (what we have been calling a gradient) and the Jacobian of these equations (what we
+ * have been calling the hessian).  We don't care if there is an underlying function that we are the gradient of.  And if we are using this
+ * to build a GLM link we are assuming the probability model is such that the gradient of log(P(y|f(obs.dot(beta)))) is our given balance equations.
+ * Note we only use the first two positions of the link.
+ * @author johnmount
+ *
+ */
+public final class BalanceBasedJacobian implements BalanceJacobianCalc {
+	public final Link link;
+	
+	public BalanceBasedJacobian(final Link link) {
+		this.link = link;
+	}
+
+	@Override
+	public BalanceJacobianCoef calc(final Obs obs, final double[] beta) {
+		final double[] z = new double[3];
+		final double bx = obs.dot(beta);
+		link.invLink(bx,z);
+		final double balanceCoef = obs.wt*(obs.y-z[0]);  // hard coded Poisson gradient component
+		final double jacobianCoef = obs.wt*(-z[1]);
+		return new BalanceJacobianCoef(balanceCoef,jacobianCoef);
+	}
+
+	@Override
+	public double evalEst(final double[] beta, final double[] x) {
+		return link.invLink(Obs.dot(beta,x));
+	}
+
+	public static final BalanceBasedJacobian poissonJacobian = new BalanceBasedJacobian(LinkBasedGradHess.LogLink);
+	public static GLMModel poissonLink = new GLMModel(poissonJacobian);
+}

Count/src/com/mzlabs/fit/BalanceJacobianCalc.java

@@ -0,0 +1,19 @@
+package com.mzlabs.fit;
+
+public interface BalanceJacobianCalc {
+	/**
+	 * 
+	 * @param obs
+	 * @param beta beta.length==obs.x.length+1
+	 * @return balance eqns and Jacobian coefficients of underlying loss function for parametes beta with respect to datum obs
+	 */
+	BalanceJacobianCoef calc(Obs obs, double[] beta);
+	
+	/**
+	 * 
+	 * @param beta
+	 * @param x x.length==beta.length-1 (last beta is the coefficient matchng an implicit constant term by convention)
+	 * @return the estimate of y given parameters beta and datum x
+	 */
+	double evalEst(double[] beta, double[] x);
+}

Count/src/com/mzlabs/fit/BalanceJacobianCoef.java

@@ -0,0 +1,20 @@
+package com.mzlabs.fit;
+
+public final class BalanceJacobianCoef {
+	public final double balanceCoef;
+	public final double jacobianCoef;
+	
+	public BalanceJacobianCoef(final double balanceCoef, final double jacobianCoef) {
+		this.balanceCoef = balanceCoef;
+		this.jacobianCoef = jacobianCoef;
+	}
+
+	public double absDiff(final BalanceJacobianCoef o) {
+		return Math.abs(balanceCoef-o.balanceCoef) + Math.abs(jacobianCoef-jacobianCoef);
+	}
+	
+	@Override
+	public String toString() {
+		return "g:" + balanceCoef + " h:" + jacobianCoef;
+	}
+}

Count/src/com/mzlabs/fit/DirectPoissonJacobian.java

@@ -0,0 +1,28 @@
+package com.mzlabs.fit;
+
+
+/**
+ * Direct implementation of gradient and Hessian coefficient calculation
+ * @author johnmount
+ *
+ */
+public final class DirectPoissonJacobian implements BalanceJacobianCalc {
+	private DirectPoissonJacobian() {
+	}
+
+	@Override
+	public BalanceJacobianCoef calc(Obs obs, double[] beta) {
+		final double bx = obs.dot(beta);
+		final double balCoef = obs.wt*(obs.y-Math.exp(bx));  // hard coded Poisson gradient component
+		final double jacCoef = obs.wt*(-Math.exp(bx));       // hard coded Poisson hessian component
+		return new BalanceJacobianCoef(balCoef,jacCoef);
+	}
+	
+	@Override
+	public double evalEst(double[] beta, double[] x) {
+		return Math.exp(Obs.dot(beta,x));
+	}
+	
+	public static final DirectPoissonJacobian poissonGradHess = new DirectPoissonJacobian();
+	public static GLMModel poissonLink = new GLMModel(poissonGradHess);
+}

Count/src/com/mzlabs/fit/GLMModel.java

@@ -4,124 +4,39 @@
 
 import com.winvector.linalg.colt.ColtMatrix;
 
-public class GLMModel implements VectorFnWithGradAndHessian {
-	public final boolean debug;
-	public final ProbYGivenZ prob;
-	public final Link link;
-	
-	
-	public GLMModel(final ProbYGivenZ prob, final Link link, final boolean debug) {
-		this.prob = prob;
-		this.link = link;
-		this.debug = debug;
-	}
-	
-	@Override
-	public double evalEst(final double[] beta, final double[] x) {
-		return link.invLink(Obs.dot(beta,x));
-	}
+public class GLMModel implements VectorFnWithJacobian {
+	public final BalanceJacobianCalc balanceJacobianCalc;
 
-	public Double checkGradCoef(final Obs obsi, final double[] beta) {
-		return null;
-	}
-	
-	public Double checkHessianCoef(final Obs obsi, final double[] beta) {
-		return null;
+	public GLMModel(final BalanceJacobianCalc balanceJacobianCalc) {
+		this.balanceJacobianCalc = balanceJacobianCalc;
 	}
 
-	public static final ProbYGivenZ PoissonProbability = new ProbYGivenZ() {
-		@Override
-		public void eval(final double y, final double z, final double[] res) {
-			if(y>0) {
-				final double fyz = Math.exp(y*Math.log(z)-z); // ignoring a gamma(y+1) term here
-				res[0] = fyz;
-				res[1] = fyz*(y/z-1);
-				res[2] = fyz*((y/z-1)*(y/z-1)-y/(z*z));
-			} else {
-				final double fyz = Math.exp(-z); // ignoring a gamma(y+1) term here
-				res[0] = fyz;
-				res[1] = -fyz;
-				res[2] = fyz;
-			}
-		}
-	};
-	public static final Link LogLink = new Link() {
-		@Override
-		public void invLink(final double z, final double[] res) {
-			final double ez = Math.exp(z);
-			res[0] = ez;
-			res[1] = ez;
-			res[2] = ez;
-		}
-		@Override
-		public double invLink(final double z) {
-			return Math.exp(z);
-		}
-	};
-	
-	public static GLMModel PoissonLink = new GLMModel(PoissonProbability,LogLink,false);
-	public static GLMModel PoissonLinkDebug = new GLMModel(PoissonProbability,LogLink,true) {
-		@Override
-		public Double checkGradCoef(final Obs obsi, final double[] beta) {
-			final double bx = obsi.dot(beta);
-			return obsi.wt*(obsi.y-Math.exp(bx));  // hard coded Poisson gradient component
-		}
-		@Override
-		public Double checkHessianCoef(final Obs obsi, final double[] beta) {
-			final double bx = obsi.dot(beta);
-			return obsi.wt*(-Math.exp(bx)); // hard coded Poisson hessian component
-		}
-	};
-
-
-	private static double sq(final double z) {
-		return z*z;
-	}
-
-	
 	@Override
-	public double lossAndGradAndHessian(final Iterable<Obs> obs, final double[] beta,
-			final double[] grad, final ColtMatrix hessian) {
+	public void balanceAndJacobian(final Iterable<Obs> obs, final double[] beta,
+			final double[] balance, final ColtMatrix jacobian) {
 		final int dim = beta.length;
-		Arrays.fill(grad,0.0);
+		Arrays.fill(balance,0.0);
 		for(int i=0;i<dim;++i) {
 			for(int j=0;j<dim;++j) {
-				hessian.set(i,j,0.0);
+				jacobian.set(i,j,0.0);
 			}
 		}
-		final double[] p = new double[3];
-		final double[] z = new double[3];
-		double sum = 0.0;
 		for(final Obs obsi: obs) {
-			final double bx = obsi.dot(beta);
-			link.invLink(bx,z);
-			prob.eval(obsi.y,z[0],p);
-			sum += obsi.wt*Math.log(z[0]);
-			final double gradScale = obsi.wt*z[1]*p[1]/p[0];
-			final double hessCoef = obsi.wt*(-sq(z[1])*sq(p[1])/sq(p[0]) + sq(z[1])*p[2]/p[0] + z[2]*p[1]/p[0]);
-			if(debug) {
-				final Double gradP = checkGradCoef(obsi,beta);
-				final Double hessP = checkHessianCoef(obsi,beta);
-				if((null!=gradP)&&(Math.abs(gradScale-gradP)>=1.0e-6)) {
-					throw new IllegalStateException("gradient checks didn't match");
-				}
-				if((null!=hessP)&&(Math.abs(hessCoef-hessP)>1.0e-6)) {
-					throw new IllegalStateException("Hessian checks didn't match");
-				}
-			}
+			final BalanceJacobianCoef ghc = balanceJacobianCalc.calc(obsi,beta);
 			for(int i=0;i<dim;++i) {
 				final double xi = i<dim-1?obsi.x[i]:1.0;
-				grad[i] += gradScale*xi;
+				balance[i] += ghc.balanceCoef*xi;
 				for(int j=0;j<dim;++j) {
 					final double xj = j<dim-1?obsi.x[j]:1.0;
-					final double hij = hessian.get(i,j);
-					hessian.set(i,j,hij+xi*xj*hessCoef);
+					final double hij = jacobian.get(i,j);
+					jacobian.set(i,j,hij+xi*xj*ghc.jacobianCoef);
 				}
 			}
 		}
-		return sum;
 	}
 
-
-
+	@Override
+	public double evalEst(double[] beta, double[] x) {
+		return balanceJacobianCalc.evalEst(beta,x);
+	}
 }

Count/src/com/mzlabs/fit/LinkBasedGradHess.java

@@ -0,0 +1,76 @@
+package com.mzlabs.fit;
+
+/**
+ * calculate the scalar portion of gradient and hessian (balance and jacobian in larger scope) directly using
+ * link and probability model
+ * @author johnmount
+ *
+ */
+public final class LinkBasedGradHess implements BalanceJacobianCalc {
+	public final ProbYGivenZ prob;
+	public final Link link;
+	
+	public LinkBasedGradHess(final ProbYGivenZ prob, final Link link) {
+		this.prob = prob;
+		this.link = link;
+	}
+
+	private static double sq(final double z) {
+		return z*z;
+	}
+	
+	@Override
+	public BalanceJacobianCoef calc(Obs obs, double[] beta) {
+		final double[] p = new double[3];
+		final double[] z = new double[3];
+		final double bx = obs.dot(beta);
+		link.invLink(bx,z);
+		prob.eval(obs.y,z[0],p);
+		final double gradCoef = obs.wt*z[1]*p[1]/p[0];
+		final double hessCoef = obs.wt*(-sq(z[1])*sq(p[1])/sq(p[0]) + sq(z[1])*p[2]/p[0] + z[2]*p[1]/p[0]);
+		return new BalanceJacobianCoef(gradCoef,hessCoef);
+	}
+	
+	@Override
+	public double evalEst(double[] beta, double[] x) {
+		return link.invLink(Obs.dot(beta,x));
+	}
+	
+	
+	
+	
+	public static final ProbYGivenZ PoissonProbability = new ProbYGivenZ() {
+		@Override
+		public void eval(final double y, final double z, final double[] res) {
+			if(y>0) {
+				final double fyz = Math.exp(y*Math.log(z)-z); // ignoring a gamma(y+1) term here
+				res[0] = fyz;
+				res[1] = fyz*(y/z-1);
+				res[2] = fyz*((y/z-1)*(y/z-1)-y/(z*z));
+			} else {
+				final double fyz = Math.exp(-z); // ignoring a gamma(y+1) term here
+				res[0] = fyz;
+				res[1] = -fyz;
+				res[2] = fyz;
+			}
+		}
+	};
+	
+	public static final Link LogLink = new Link() {
+		@Override
+		public void invLink(final double z, final double[] res) {
+			final double ez = Math.exp(z);
+			res[0] = ez;
+			res[1] = ez;
+			res[2] = ez;
+		}
+		@Override
+		public double invLink(final double z) {
+			return Math.exp(z);
+		}
+	};
+	
+	
+	public static final LinkBasedGradHess poissonGradHess = new LinkBasedGradHess(PoissonProbability,LogLink);
+	public static GLMModel poissonLink = new GLMModel(poissonGradHess);
+}