detect numerical problems in Q.R decomposition for Levenberg-Marquardt estimator

and report them appropriately
JIRA: MATH-199

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@640204 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-03-23 13:36:03 +00:00
parent cd5f65c4c7
commit faead6c3fb
3 changed files with 92 additions and 1 deletions

View File

@ -733,8 +733,9 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
* are performed in non-increasing columns norms order thanks to columns
* pivoting. The diagonal elements of the R matrix are therefore also in
* non-increasing absolute values order.</p>
* @exception EstimationException if the decomposition cannot be performed
*/
private void qrDecomposition() {
private void qrDecomposition() throws EstimationException {
// initializations
for (int k = 0; k < cols; ++k) {
@ -760,6 +761,10 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
double aki = jacobian[index];
norm2 += aki * aki;
}
if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
throw new EstimationException("unable to perform Q.R decomposition on the {0}x{1} jacobian matrix",
new Object[] { new Integer(rows), new Integer(cols) });
}
if (norm2 > ak2) {
nextColumn = i;
ak2 = norm2;

View File

@ -47,6 +47,10 @@ Commons Math Release Notes</title>
<action dev="luc" type="fix" issue="MATH-198" due-to="Frederick Salardi">
added an error detection for missing imaginary character while parsing complex string
</action>
<action dev="luc" type="fix" issue="MATH-199" due-to="Mick">
detect numerical problems in Q.R decomposition for Levenberg-Marquardt estimator
and report them appropriately
</action>
</release>
<release version="1.2" date="2008-02-24"
description="This release combines bug fixes and new features. Most notable

View File

@ -20,6 +20,7 @@ package org.apache.commons.math.estimation;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.Iterator;
import java.util.Locale;
import java.util.Set;
import org.apache.commons.math.estimation.EstimatedParameter;
@ -587,6 +588,22 @@ public class LevenbergMarquardtEstimatorTest
assertEquals( 0.20750021499570379, circle.getY(), 1.0e-8);
}
public void testMath199() {
try {
QuadraticProblem problem = new QuadraticProblem();
problem.addPoint (0, -3.182591015485607, 0.0);
problem.addPoint (1, -2.5581184967730577, 4.4E-323);
problem.addPoint (2, -2.1488478161387325, 1.0);
problem.addPoint (3, -1.9122489313410047, 4.4E-323);
problem.addPoint (4, 1.7785661310051026, 0.0);
new LevenbergMarquardtEstimator().estimate(problem);
fail("an exception should have been thrown");
} catch (EstimationException ee) {
// expected behavior
}
}
private static class LinearProblem implements EstimationProblem {
public LinearProblem(LinearMeasurement[] measurements) {
@ -759,6 +776,71 @@ public class LevenbergMarquardtEstimatorTest
private ArrayList points;
}
public class QuadraticProblem extends SimpleEstimationProblem {
private EstimatedParameter a;
private EstimatedParameter b;
private EstimatedParameter c;
public QuadraticProblem() {
a = new EstimatedParameter("a", 0.0);
b = new EstimatedParameter("b", 0.0);
c = new EstimatedParameter("c", 0.0);
addParameter(a);
addParameter(b);
addParameter(c);
}
public void addPoint(double x, double y, double w) {
addMeasurement(new LocalMeasurement(x, y, w));
}
public double getA() {
return a.getEstimate();
}
public double getB() {
return b.getEstimate();
}
public double getC() {
return c.getEstimate();
}
public double theoreticalValue(double x) {
return ( (a.getEstimate() * x + b.getEstimate() ) * x + c.getEstimate());
}
private double partial(double x, EstimatedParameter parameter) {
if (parameter == a) {
return x * x;
} else if (parameter == b) {
return x;
} else {
return 1.0;
}
}
private class LocalMeasurement extends WeightedMeasurement {
private final double x;
// constructor
public LocalMeasurement(double x, double y, double w) {
super(w, y);
this.x = x;
}
public double getTheoreticalValue() {
return theoreticalValue(x);
}
public double getPartial(EstimatedParameter parameter) {
return partial(x, parameter);
}
}
}
public static Test suite() {
return new TestSuite(LevenbergMarquardtEstimatorTest.class);