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");
}