/**
* Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
*
* Please see distribution for license.
*/
package com.opengamma.analytics.math.statistics.leastsquare;
import java.util.Arrays;
import org.apache.commons.lang.Validate;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import com.opengamma.analytics.math.FunctionUtils;
import com.opengamma.analytics.math.MathException;
import com.opengamma.analytics.math.differentiation.VectorFieldFirstOrderDifferentiator;
import com.opengamma.analytics.math.differentiation.VectorFieldSecondOrderDifferentiator;
import com.opengamma.analytics.math.function.Function1D;
import com.opengamma.analytics.math.function.ParameterizedFunction;
import com.opengamma.analytics.math.linearalgebra.Decomposition;
import com.opengamma.analytics.math.linearalgebra.DecompositionFactory;
import com.opengamma.analytics.math.linearalgebra.DecompositionResult;
import com.opengamma.analytics.math.linearalgebra.SVDecompositionCommons;
import com.opengamma.analytics.math.linearalgebra.SVDecompositionResult;
import com.opengamma.analytics.math.matrix.DoubleMatrix1D;
import com.opengamma.analytics.math.matrix.DoubleMatrix2D;
import com.opengamma.analytics.math.matrix.DoubleMatrixUtils;
import com.opengamma.analytics.math.matrix.MatrixAlgebra;
import com.opengamma.analytics.math.matrix.MatrixAlgebraFactory;
import com.opengamma.util.ArgumentChecker;
/**
*
*/
public class NonLinearLeastSquare {
private static final Logger LOGGER = LoggerFactory.getLogger(NonLinearLeastSquare.class);
private static final int MAX_ATTEMPTS = 10000;
private static final Function1D<DoubleMatrix1D, Boolean> UNCONSTAINED = new Function1D<DoubleMatrix1D, Boolean>() {
@Override
public Boolean evaluate(final DoubleMatrix1D x) {
return true;
}
};
private final double _eps;
private final Decomposition<?> _decomposition;
private final MatrixAlgebra _algebra;
public NonLinearLeastSquare() {
this(DecompositionFactory.SV_COMMONS, MatrixAlgebraFactory.OG_ALGEBRA, 1e-8);
}
public NonLinearLeastSquare(final Decomposition<?> decomposition, final MatrixAlgebra algebra, final double eps) {
_decomposition = decomposition;
_algebra = algebra;
_eps = eps;
}
/**
* Use this when the model is in the ParameterizedFunction form and analytic parameter sensitivity is not available
* @param x Set of measurement points
* @param y Set of measurement values
* @param func The model in ParameterizedFunction form (i.e. takes measurement points and a set of parameters and returns a model value)
* @param startPos Initial value of the parameters
* @return A LeastSquareResults object
*/
public LeastSquareResults solve(final DoubleMatrix1D x, final DoubleMatrix1D y, final ParameterizedFunction<Double, DoubleMatrix1D, Double> func,
final DoubleMatrix1D startPos) {
Validate.notNull(x, "x");
Validate.notNull(y, "y");
final int n = x.getNumberOfElements();
Validate.isTrue(y.getNumberOfElements() == n, "y wrong length");
final double[] sigmas = new double[n];
Arrays.fill(sigmas, 1); //emcleod 31-1-2011 arbitrary value for now
return solve(x, y, new DoubleMatrix1D(sigmas), func, startPos);
}
/**
* Use this when the model is in the ParameterizedFunction form and analytic parameter sensitivity is not available but a measurement error is.
* @param x Set of measurement points
* @param y Set of measurement values
* @param sigma y Set of measurement errors
* @param func The model in ParameterizedFunction form (i.e. takes measurement points and a set of parameters and returns a model value)
* @param startPos Initial value of the parameters
* @return A LeastSquareResults object
*/
public LeastSquareResults solve(final DoubleMatrix1D x, final DoubleMatrix1D y, final double sigma, final ParameterizedFunction<Double, DoubleMatrix1D, Double> func,
final DoubleMatrix1D startPos) {
Validate.notNull(x, "x");
Validate.notNull(y, "y");
Validate.notNull(sigma, "sigma");
final int n = x.getNumberOfElements();
Validate.isTrue(y.getNumberOfElements() == n, "y wrong length");
final double[] sigmas = new double[n];
Arrays.fill(sigmas, sigma);
return solve(x, y, new DoubleMatrix1D(sigmas), func, startPos);
}
/**
* Use this when the model is in the ParameterizedFunction form and analytic parameter sensitivity is not available but an array of measurements errors is.
* @param x Set of measurement points
* @param y Set of measurement values
* @param sigma Set of measurement errors
* @param func The model in ParameterizedFunction form (i.e. takes measurement points and a set of parameters and returns a model value)
* @param startPos Initial value of the parameters
* @return A LeastSquareResults object
*/
public LeastSquareResults solve(final DoubleMatrix1D x, final DoubleMatrix1D y, final DoubleMatrix1D sigma,
final ParameterizedFunction<Double, DoubleMatrix1D, Double> func, final DoubleMatrix1D startPos) {
Validate.notNull(x, "x");
Validate.notNull(y, "y");
Validate.notNull(sigma, "sigma");
final int n = x.getNumberOfElements();
Validate.isTrue(y.getNumberOfElements() == n, "y wrong length");
Validate.isTrue(sigma.getNumberOfElements() == n, "sigma wrong length");
final Function1D<DoubleMatrix1D, DoubleMatrix1D> func1D = new Function1D<DoubleMatrix1D, DoubleMatrix1D>() {
@Override
public DoubleMatrix1D evaluate(final DoubleMatrix1D theta) {
final int m = x.getNumberOfElements();
final double[] res = new double[m];
for (int i = 0; i < m; i++) {
res[i] = func.evaluate(x.getEntry(i), theta);
}
return new DoubleMatrix1D(res);
}
};
return solve(y, sigma, func1D, startPos, null);
}
/**
* Use this when the model is in the ParameterizedFunction form and analytic parameter sensitivity
* @param x Set of measurement points
* @param y Set of measurement values
* @param func The model in ParameterizedFunction form (i.e. takes a measurement points and a set of parameters and returns a model value)
* @param grad The model parameter sensitivities in ParameterizedFunction form (i.e. takes a measurement points and a set of parameters and returns a model parameter sensitivities)
* @param startPos Initial value of the parameters
* @return value of the fitted parameters
*/
public LeastSquareResults solve(final DoubleMatrix1D x, final DoubleMatrix1D y, final ParameterizedFunction<Double, DoubleMatrix1D, Double> func,
final ParameterizedFunction<Double, DoubleMatrix1D, DoubleMatrix1D> grad, final DoubleMatrix1D startPos) {
Validate.notNull(x, "x");
Validate.notNull(y, "y");
Validate.notNull(x, "sigma");
final int n = x.getNumberOfElements();
Validate.isTrue(y.getNumberOfElements() == n, "y wrong length");
final double[] sigmas = new double[n];
Arrays.fill(sigmas, 1); //emcleod 31-1-2011 arbitrary value for now
return solve(x, y, new DoubleMatrix1D(sigmas), func, grad, startPos);
}
/**
* Use this when the model is in the ParameterizedFunction form and analytic parameter sensitivity and a single measurement error are available
* @param x Set of measurement points
* @param y Set of measurement values
* @param sigma Measurement errors
* @param func The model in ParameterizedFunction form (i.e. takes a measurement points and a set of parameters and returns a model value)
* @param grad The model parameter sensitivities in ParameterizedFunction form (i.e. takes a measurement points and a set of parameters and returns a model parameter sensitivities)
* @param startPos Initial value of the parameters
* @return value of the fitted parameters
*/
public LeastSquareResults solve(final DoubleMatrix1D x, final DoubleMatrix1D y, final double sigma, final ParameterizedFunction<Double, DoubleMatrix1D, Double> func,
final ParameterizedFunction<Double, DoubleMatrix1D, DoubleMatrix1D> grad, final DoubleMatrix1D startPos) {
Validate.notNull(x, "x");
Validate.notNull(y, "y");
final int n = x.getNumberOfElements();
Validate.isTrue(y.getNumberOfElements() == n, "y wrong length");
final double[] sigmas = new double[n];
Arrays.fill(sigmas, sigma);
return solve(x, y, new DoubleMatrix1D(sigmas), func, grad, startPos);
}
/**
* Use this when the model is in the ParameterizedFunction form and analytic parameter sensitivity and measurement errors are available
* @param x Set of measurement points
* @param y Set of measurement values
* @param sigma Set of measurement errors
* @param func The model in ParameterizedFunction form (i.e. takes a measurement points and a set of parameters and returns a model value)
* @param grad The model parameter sensitivities in ParameterizedFunction form (i.e. takes a measurement points and a set of parameters and returns a model parameter sensitivities)
* @param startPos Initial value of the parameters
* @return value of the fitted parameters
*/
public LeastSquareResults solve(final DoubleMatrix1D x, final DoubleMatrix1D y, final DoubleMatrix1D sigma,
final ParameterizedFunction<Double, DoubleMatrix1D, Double> func, final ParameterizedFunction<Double, DoubleMatrix1D, DoubleMatrix1D> grad,
final DoubleMatrix1D startPos) {
Validate.notNull(x, "x");
Validate.notNull(y, "y");
Validate.notNull(x, "sigma");
final int n = x.getNumberOfElements();
Validate.isTrue(y.getNumberOfElements() == n, "y wrong length");
Validate.isTrue(sigma.getNumberOfElements() == n, "sigma wrong length");
final Function1D<DoubleMatrix1D, DoubleMatrix1D> func1D = new Function1D<DoubleMatrix1D, DoubleMatrix1D>() {
@Override
public DoubleMatrix1D evaluate(final DoubleMatrix1D theta) {
final int m = x.getNumberOfElements();
final double[] res = new double[m];
for (int i = 0; i < m; i++) {
res[i] = func.evaluate(x.getEntry(i), theta);
}
return new DoubleMatrix1D(res);
}
};
final Function1D<DoubleMatrix1D, DoubleMatrix2D> jac = new Function1D<DoubleMatrix1D, DoubleMatrix2D>() {
@Override
public DoubleMatrix2D evaluate(final DoubleMatrix1D theta) {
final int m = x.getNumberOfElements();
final double[][] res = new double[m][];
for (int i = 0; i < m; i++) {
final DoubleMatrix1D temp = grad.evaluate(x.getEntry(i), theta);
res[i] = temp.getData();
}
return new DoubleMatrix2D(res);
}
};
return solve(y, sigma, func1D, jac, startPos, null);
}
/**
* Use this when the model is given as a function of its parameters only (i.e. a function that takes a set of parameters and return a set of model values,
* so the measurement points are already known to the function), and analytic parameter sensitivity is not available
* @param observedValues Set of measurement values
* @param func The model as a function of its parameters only
* @param startPos Initial value of the parameters
* @return value of the fitted parameters
*/
public LeastSquareResults solve(final DoubleMatrix1D observedValues, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func, final DoubleMatrix1D startPos) {
final int n = observedValues.getNumberOfElements();
final VectorFieldFirstOrderDifferentiator jac = new VectorFieldFirstOrderDifferentiator();
return solve(observedValues, new DoubleMatrix1D(n, 1.0), func, jac.differentiate(func), startPos, null);
}
/**
* Use this when the model is given as a function of its parameters only (i.e. a function that takes a set of parameters and return a set of model values,
* so the measurement points are already known to the function), and analytic parameter sensitivity is not available
* @param observedValues Set of measurement values
* @param sigma Set of measurement errors
* @param func The model as a function of its parameters only
* @param startPos Initial value of the parameters
* @return value of the fitted parameters
*/
public LeastSquareResults solve(final DoubleMatrix1D observedValues, final DoubleMatrix1D sigma, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func,
final DoubleMatrix1D startPos) {
final VectorFieldFirstOrderDifferentiator jac = new VectorFieldFirstOrderDifferentiator();
return solve(observedValues, sigma, func, jac.differentiate(func), startPos, null);
}
/**
* Use this when the model is given as a function of its parameters only (i.e. a function that takes a set of parameters and return a set of model values,
* so the measurement points are already known to the function), and analytic parameter sensitivity is not available
* @param observedValues Set of measurement values
* @param sigma Set of measurement errors
* @param func The model as a function of its parameters only
* @param startPos Initial value of the parameters
* @param maxJumps A vector containing the maximum absolute allowed step in a particular direction in each iteration. Can be null, in which case no constant
* on the step size is applied.
* @return value of the fitted parameters
*/
public LeastSquareResults solve(final DoubleMatrix1D observedValues, final DoubleMatrix1D sigma, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func,
final DoubleMatrix1D startPos, final DoubleMatrix1D maxJumps) {
final VectorFieldFirstOrderDifferentiator jac = new VectorFieldFirstOrderDifferentiator();
return solve(observedValues, sigma, func, jac.differentiate(func), startPos, maxJumps);
}
/**
* Use this when the model is given as a function of its parameters only (i.e. a function that takes a set of parameters and return a set of model values,
* so the measurement points are already known to the function), and analytic parameter sensitivity is available
* @param observedValues Set of measurement values
* @param sigma Set of measurement errors
* @param func The model as a function of its parameters only
* @param jac The model sensitivity to its parameters (i.e. the Jacobian matrix) as a function of its parameters only
* @param startPos Initial value of the parameters
* @return value of the fitted parameters
*/
public LeastSquareResults solve(final DoubleMatrix1D observedValues, final DoubleMatrix1D sigma, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func,
final Function1D<DoubleMatrix1D, DoubleMatrix2D> jac, final DoubleMatrix1D startPos) {
return solve(observedValues, sigma, func, jac, startPos, UNCONSTAINED, null);
}
/**
* Use this when the model is given as a function of its parameters only (i.e. a function that takes a set of parameters and return a set of model values,
* so the measurement points are already known to the function), and analytic parameter sensitivity is available
* @param observedValues Set of measurement values
* @param sigma Set of measurement errors
* @param func The model as a function of its parameters only
* @param jac The model sensitivity to its parameters (i.e. the Jacobian matrix) as a function of its parameters only
* @param startPos Initial value of the parameters
* @param maxJumps A vector containing the maximum absolute allowed step in a particular direction in each iteration. Can be null, in which case on constant
* on the step size is applied.
* @return value of the fitted parameters
*/
public LeastSquareResults solve(final DoubleMatrix1D observedValues, final DoubleMatrix1D sigma, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func,
final Function1D<DoubleMatrix1D, DoubleMatrix2D> jac, final DoubleMatrix1D startPos, final DoubleMatrix1D maxJumps) {
return solve(observedValues, sigma, func, jac, startPos, UNCONSTAINED, maxJumps);
}
/**
* Use this when the model is given as a function of its parameters only (i.e. a function that takes a set of parameters and return a set of model values,
* so the measurement points are already known to the function), and analytic parameter sensitivity is available
* @param observedValues Set of measurement values
* @param sigma Set of measurement errors
* @param func The model as a function of its parameters only
* @param jac The model sensitivity to its parameters (i.e. the Jacobian matrix) as a function of its parameters only
* @param startPos Initial value of the parameters
* @param constraints A function that returns true if the trial point is within the constraints of the model
* @param maxJumps A vector containing the maximum absolute allowed step in a particular direction in each iteration. Can be null, in which case on constant
* on the step size is applied.
* @return value of the fitted parameters
*/
public LeastSquareResults solve(final DoubleMatrix1D observedValues, final DoubleMatrix1D sigma, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func,
final Function1D<DoubleMatrix1D, DoubleMatrix2D> jac, final DoubleMatrix1D startPos, final Function1D<DoubleMatrix1D, Boolean> constraints,
final DoubleMatrix1D maxJumps) {
Validate.notNull(observedValues, "observedValues");
Validate.notNull(sigma, " sigma");
Validate.notNull(func, " func");
Validate.notNull(jac, " jac");
Validate.notNull(startPos, "startPos");
final int nObs = observedValues.getNumberOfElements();
final int nParms = startPos.getNumberOfElements();
Validate.isTrue(nObs == sigma.getNumberOfElements(), "observedValues and sigma must be same length");
ArgumentChecker.isTrue(nObs >= nParms, "must have data points greater or equal to number of parameters. #date points = {}, #parameters = {}", nObs, nParms);
ArgumentChecker.isTrue(constraints.evaluate(startPos), "The inital value of the parameters (startPos) is {} - this is not an allowed value", startPos);
DoubleMatrix2D alpha;
DecompositionResult decmp;
DoubleMatrix1D theta = startPos;
double lambda = 0.0; //TODO debug if the model is linear, it will be solved in 1 step
double newChiSqr, oldChiSqr;
DoubleMatrix1D error = getError(func, observedValues, sigma, theta);
DoubleMatrix1D newError;
DoubleMatrix2D jacobian = getJacobian(jac, sigma, theta);
oldChiSqr = getChiSqr(error);
//If we start at the solution we are done
if (oldChiSqr == 0.0) {
return finish(oldChiSqr, jacobian, theta, sigma);
}
DoubleMatrix1D beta = getChiSqrGrad(error, jacobian);
for (int count = 0; count < MAX_ATTEMPTS; count++) {
alpha = getModifiedCurvatureMatrix(jacobian, lambda);
DoubleMatrix1D deltaTheta;
try {
decmp = _decomposition.evaluate(alpha);
deltaTheta = decmp.solve(beta);
} catch (final Exception e) {
throw new MathException(e);
}
DoubleMatrix1D trialTheta = (DoubleMatrix1D) _algebra.add(theta, deltaTheta);
//if the new value of theta is not in the model domain or the jump is too large, keep increasing lambda until an acceptable step is found
if (!constraints.evaluate(trialTheta) || !allowJump(deltaTheta, maxJumps)) {
lambda = increaseLambda(lambda);
continue;
}
newError = getError(func, observedValues, sigma, trialTheta);
newChiSqr = getChiSqr(newError);
//Check for convergence when no improvement in chiSqr occurs
if (Math.abs(newChiSqr - oldChiSqr) / (1 + oldChiSqr) < _eps) {
final DoubleMatrix2D alpha0 = lambda == 0.0 ? alpha : getModifiedCurvatureMatrix(jacobian, 0.0);
//if the model is an exact fit to the data, then no more improvement is possible
if (newChiSqr < _eps) {
if (lambda > 0.0) {
decmp = _decomposition.evaluate(alpha0);
}
return finish(alpha0, decmp, newChiSqr, jacobian, trialTheta, sigma);
}
final SVDecompositionCommons svd = (SVDecompositionCommons) DecompositionFactory.SV_COMMONS;
//add the second derivative information to the Hessian matrix to check we are not at a local maximum or saddle point
final VectorFieldSecondOrderDifferentiator diff = new VectorFieldSecondOrderDifferentiator();
final Function1D<DoubleMatrix1D, DoubleMatrix2D[]> secDivFunc = diff.differentiate(func, constraints);
final DoubleMatrix2D[] secDiv = secDivFunc.evaluate(trialTheta);
final double[][] temp = new double[nParms][nParms];
for (int i = 0; i < nObs; i++) {
for (int j = 0; j < nParms; j++) {
for (int k = 0; k < nParms; k++) {
temp[j][k] -= newError.getEntry(i) * secDiv[i].getEntry(j, k) / sigma.getEntry(i);
}
}
}
final DoubleMatrix2D newAlpha = (DoubleMatrix2D) _algebra.add(alpha0, new DoubleMatrix2D(temp));
final SVDecompositionResult svdRes = svd.evaluate(newAlpha);
final double[] w = svdRes.getSingularValues();
final DoubleMatrix2D u = svdRes.getU();
final DoubleMatrix2D v = svdRes.getV();
final double[] p = new double[nParms];
boolean saddle = false;
double sum = 0.0;
for (int i = 0; i < nParms; i++) {
double a = 0.0;
for (int j = 0; j < nParms; j++) {
a += u.getEntry(j, i) * v.getEntry(j, i);
}
final int sign = a > 0.0 ? 1 : -1;
if (w[i] * sign < 0.0) {
sum += w[i];
w[i] = -w[i];
saddle = true;
}
}
//if a local maximum or saddle point is found (as indicated by negative eigenvalues), move in a direction that is a weighted
//sum of the eigenvectors corresponding to the negative eigenvalues
if (saddle) {
lambda = increaseLambda(lambda);
for (int i = 0; i < nParms; i++) {
if (w[i] < 0.0) {
final double scale = 0.5 * Math.sqrt(-oldChiSqr * w[i]) / sum;
for (int j = 0; j < nParms; j++) {
p[j] += scale * u.getEntry(j, i);
}
}
}
final DoubleMatrix1D direction = new DoubleMatrix1D(p);
deltaTheta = direction;
trialTheta = (DoubleMatrix1D) _algebra.add(theta, deltaTheta);
int i = 0;
double scale = 1.0;
while (!constraints.evaluate(trialTheta)) {
scale *= -0.5;
deltaTheta = (DoubleMatrix1D) _algebra.scale(direction, scale);
trialTheta = (DoubleMatrix1D) _algebra.add(theta, deltaTheta);
i++;
if (i > 10) {
throw new MathException("Could not satify constraint");
}
}
newError = getError(func, observedValues, sigma, trialTheta);
newChiSqr = getChiSqr(newError);
int counter = 0;
while (newChiSqr > oldChiSqr) {
//if even a tiny move along the negative eigenvalue cannot improve chiSqr, then exit
if (counter > 10 || Math.abs(newChiSqr - oldChiSqr) / (1 + oldChiSqr) < _eps) {
LOGGER.warn("Saddle point detected, but no improvement to chi^2 possible by moving away. It is recommended that a different starting point is used.");
return finish(newAlpha, decmp, oldChiSqr, jacobian, theta, sigma);
}
scale /= 2.0;
deltaTheta = (DoubleMatrix1D) _algebra.scale(direction, scale);
trialTheta = (DoubleMatrix1D) _algebra.add(theta, deltaTheta);
newError = getError(func, observedValues, sigma, trialTheta);
newChiSqr = getChiSqr(newError);
counter++;
}
} else {
//this should be the normal finish - i.e. no improvement in chiSqr and at a true minimum (although there is no guarantee it is not a local minimum)
return finish(newAlpha, decmp, newChiSqr, jacobian, trialTheta, sigma);
}
}
if (newChiSqr < oldChiSqr) {
lambda = decreaseLambda(lambda);
theta = trialTheta;
error = newError;
jacobian = getJacobian(jac, sigma, trialTheta);
beta = getChiSqrGrad(error, jacobian);
// check for convergence
// if (_algebra.getNorm2(beta) < _eps * g0) {
// return finish(newChiSqr, jacobian, trialTheta, sigma);
// }
oldChiSqr = newChiSqr;
} else {
lambda = increaseLambda(lambda);
}
}
throw new MathException("Could not converge in " + MAX_ATTEMPTS + " attempts");
}
private double decreaseLambda(final double lambda) {
return lambda / 10;
}
private double increaseLambda(final double lambda) {
if (lambda == 0.0) { // this will happen the first time a full quadratic step fails
return 0.1;
}
return lambda * 10;
}
private boolean allowJump(final DoubleMatrix1D deltaTheta, final DoubleMatrix1D maxJumps) {
if (maxJumps == null) {
return true;
}
final int n = deltaTheta.getNumberOfElements();
for (int i = 0; i < n; i++) {
if (Math.abs(deltaTheta.getEntry(i)) > maxJumps.getEntry(i)) {
return false;
}
}
return true;
}
/**
*
* the inverse-Jacobian where the i-j entry is the sensitivity of the ith (fitted) parameter (a_i) to the jth data point (y_j).
* @param sigma Set of measurement errors
* @param func The model as a function of its parameters only
* @param jac The model sensitivity to its parameters (i.e. the Jacobian matrix) as a function of its parameters only
* @param originalSolution The value of the parameters at a converged solution
* @return inverse-Jacobian
*/
public DoubleMatrix2D calInverseJacobian(final DoubleMatrix1D sigma, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func,
final Function1D<DoubleMatrix1D, DoubleMatrix2D> jac, final DoubleMatrix1D originalSolution) {
final DoubleMatrix2D jacobian = getJacobian(jac, sigma, originalSolution);
final DoubleMatrix2D a = getModifiedCurvatureMatrix(jacobian, 0.0);
final DoubleMatrix2D bT = getBTranspose(jacobian, sigma);
final DecompositionResult decRes = _decomposition.evaluate(a);
return decRes.solve(bT);
}
private LeastSquareResults finish(final double newChiSqr, final DoubleMatrix2D jacobian, final DoubleMatrix1D newTheta, final DoubleMatrix1D sigma) {
final DoubleMatrix2D alpha = getModifiedCurvatureMatrix(jacobian, 0.0);
final DecompositionResult decmp = _decomposition.evaluate(alpha);
return finish(alpha, decmp, newChiSqr, jacobian, newTheta, sigma);
}
private LeastSquareResults finish(final DoubleMatrix2D alpha, final DecompositionResult decmp, final double newChiSqr, final DoubleMatrix2D jacobian,
final DoubleMatrix1D newTheta, final DoubleMatrix1D sigma) {
final DoubleMatrix2D covariance = decmp.solve(DoubleMatrixUtils.getIdentityMatrix2D(alpha.getNumberOfRows()));
final DoubleMatrix2D bT = getBTranspose(jacobian, sigma);
final DoubleMatrix2D inverseJacobian = decmp.solve(bT);
return new LeastSquareResults(newChiSqr, newTheta, covariance, inverseJacobian);
}
private DoubleMatrix1D getError(final Function1D<DoubleMatrix1D, DoubleMatrix1D> func, final DoubleMatrix1D observedValues, final DoubleMatrix1D sigma,
final DoubleMatrix1D theta) {
final int n = observedValues.getNumberOfElements();
final DoubleMatrix1D modelValues = func.evaluate(theta);
Validate.isTrue(n == modelValues.getNumberOfElements(), "Number of data points different between model (" + modelValues.getNumberOfElements() + ") and observed ("
+ n + ")");
final double[] res = new double[n];
for (int i = 0; i < n; i++) {
res[i] = (observedValues.getEntry(i) - modelValues.getEntry(i)) / sigma.getEntry(i);
}
return new DoubleMatrix1D(res);
}
private DoubleMatrix2D getBTranspose(final DoubleMatrix2D jacobian, final DoubleMatrix1D sigma) {
final int n = jacobian.getNumberOfRows();
final int m = jacobian.getNumberOfColumns();
final double[][] res = new double[m][n];
for (int k = 0; k < m; k++) {
for (int i = 0; i < n; i++) {
res[k][i] = jacobian.getEntry(i, k) / sigma.getEntry(i);
}
}
return new DoubleMatrix2D(res);
}
private DoubleMatrix2D getJacobian(final Function1D<DoubleMatrix1D, DoubleMatrix2D> jac, final DoubleMatrix1D sigma, final DoubleMatrix1D theta) {
final DoubleMatrix2D res = jac.evaluate(theta);
final double[][] data = res.getData();
final int n = res.getNumberOfRows();
final int m = res.getNumberOfColumns();
Validate.isTrue(theta.getNumberOfElements() == m, "Jacobian is wrong size");
Validate.isTrue(sigma.getNumberOfElements() == n, "Jacobian is wrong size");
for (int i = 0; i < n; i++) {
for (int j = 0; j < m; j++) {
data[i][j] /= sigma.getEntry(i);
}
}
return res;
}
private double getChiSqr(final DoubleMatrix1D error) {
return _algebra.getInnerProduct(error, error);
}
private DoubleMatrix1D getChiSqrGrad(final DoubleMatrix1D error, final DoubleMatrix2D jacobian) {
return (DoubleMatrix1D) _algebra.multiply(error, jacobian);
}
@SuppressWarnings("unused")
private DoubleMatrix1D getDiagonalCurvatureMatrix(final DoubleMatrix2D jacobian) {
final int n = jacobian.getNumberOfRows();
final int m = jacobian.getNumberOfColumns();
final double[] alpha = new double[m];
for (int i = 0; i < m; i++) {
double sum = 0.0;
for (int k = 0; k < n; k++) {
sum += FunctionUtils.square(jacobian.getEntry(k, i));
}
alpha[i] = sum;
}
return new DoubleMatrix1D(alpha);
}
private DoubleMatrix2D getModifiedCurvatureMatrix(final DoubleMatrix2D jacobian, final double lambda) {
final int n = jacobian.getNumberOfRows();
final int m = jacobian.getNumberOfColumns();
final double[][] alpha = new double[m][m];
for (int i = 0; i < m; i++) {
double sum = 0.0;
for (int k = 0; k < n; k++) {
sum += FunctionUtils.square(jacobian.getEntry(k, i));
}
alpha[i][i] = (1 + lambda) * sum;
for (int j = i + 1; j < m; j++) {
sum = 0.0;
for (int k = 0; k < n; k++) {
sum += jacobian.getEntry(k, i) * jacobian.getEntry(k, j);
}
alpha[i][j] = sum;
alpha[j][i] = sum;
}
}
return new DoubleMatrix2D(alpha);
}
}