commons-issues mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From "Luc Maisonobe (JIRA)" <j...@apache.org>
Subject [jira] Commented: (MATH-199) exception in LevenbergMarquardtEstimator
Date Mon, 24 Mar 2008 15:09:24 GMT

    [ https://issues.apache.org/jira/browse/MATH-199?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=12581558#action_12581558
] 

Luc Maisonobe commented on MATH-199:
------------------------------------

Yes, this is this small value that triggers the error. It is a valid number but during the
Q.R decomposition, this number is used in several operations. One of these operations lead
to an overflow.

At line 785 of LevenbergMarquardt.java, we compute :  double betak = 1.0 / (ak2 - akk * alpha)

With the very small values of the example, the value for betak overflows while processing
the second column (first column is computed without problem). The denominator ak2 - akk *
alpha has a value of 1.43e-322 which can be handled by double. It's inverse is 6.9e321, far
larger than the Double.MAX_VALUE which is about 1.8e+308. The result of the computation is
that betak is set to positive infinity. The rest of the computation behaves badly with such
a value. Some of the elements of third column are set to infinity, others are set to NaN.

The fact is that despite many representable double number have a representable inverse in
IEE754, it is not true for all.

> exception in LevenbergMarquardtEstimator
> ----------------------------------------
>
>                 Key: MATH-199
>                 URL: https://issues.apache.org/jira/browse/MATH-199
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 1.2
>         Environment: Windows XP
> Java 6
>            Reporter: Mick
>
> I get this exception:
> Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: -1
>        at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.qrDecomposition(LevenbergMarquardtEstimator.java:772)
>        at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.estimate(LevenbergMarquardtEstimator.java:232)
>        at quadraticFitterProblem.QuadraticFitterProblem.<init>(QuadraticFitterProblem.java:27)
>        at quadraticFitterProblem.QuadraticFitterProblem.main(QuadraticFitterProblem.java:40)
> on the code below.
> The exception does not occur all the weights in the quadraticFitter are 0.0;
> ---------------------------------------------------------------------------------------------
> package quadraticFitterProblem;
> import org.apache.commons.math.estimation.EstimationException;
> import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
> //import org.apache.commons.math.estimation.WeightedMeasurement;
> import com.strategicanalytics.dtd.data.smoothers.QuadraticFitter;
> public class QuadraticFitterProblem {
>        private QuadraticFitter quadraticFitter;
>        public QuadraticFitterProblem() {
>          // create the uninitialized fitting problem
>          quadraticFitter = new QuadraticFitter();
>          quadraticFitter.addPoint (0,  -3.182591015485607, 0.0);
>          quadraticFitter.addPoint (1,  -2.5581184967730577, 4.4E-323);
>          quadraticFitter.addPoint (2,  -2.1488478161387325, 1.0);
>          quadraticFitter.addPoint (3,  -1.9122489313410047, 4.4E-323);
>          quadraticFitter.addPoint (4,  1.7785661310051026, 0.0);
>          try {
>            // solve the problem, using a Levenberg-Marquardt algorithm with
> default settings
>            LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
>            //WeightedMeasurement[] wm = quadraticFitter.getMeasurements();
>            estimator.estimate(quadraticFitter);
>          } catch (EstimationException ee) {
>                System.err.println(ee.getMessage());
>          }
>        }
>        /**
>         * @param args
>         *
>         */
>        public static void main(String[] args) {
>                        new QuadraticFitterProblem();
>                        System.out.println ("Done.");
>        }
> }
> ----------------------------------------------------------------------------------------------
> import org.apache.commons.math.estimation.EstimatedParameter;
> //import org.apache.commons.math.estimation.EstimationException;
> //import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
> import org.apache.commons.math.estimation.SimpleEstimationProblem;
> import org.apache.commons.math.estimation.WeightedMeasurement;
> public class QuadraticFitter extends SimpleEstimationProblem {
>        // y = a x<sup>2</sup> + b x + c
>    private EstimatedParameter a;
>    private EstimatedParameter b;
>    private EstimatedParameter c;
>    /**
>     * constructor
>     *
>     *Fitter for a quadratic model to a sample of 2D points.
>     * <p>The model is y(x) = a x<sup>2</sup> + b x + c
>     * its three parameters of the model are a, b and c.</p>
>     */
>    public QuadraticFitter() {
>        // three parameters of the model
>        a = new EstimatedParameter("a", 0.0);
>        b = new EstimatedParameter("b", 0.0);
>        c = new EstimatedParameter("c", 0.0);
>        // provide the parameters to the base class which
>        // implements the getAllParameters and getUnboundParameters methods
>        addParameter(a);
>        addParameter(b);
>        addParameter(c);
>    }
>    /**
>     * Add a sample point
>     *
>     * @param x abscissa
>     * @param y ordinate
>     * @param w weight
>     */
>    public void addPoint(double x, double y, double w) {
>        addMeasurement(new LocalMeasurement(x, y, w));
>    }
>    /**
>     * Get the value of the quadratic coefficient.
>     *
>     * @return the value of a for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getA() {
>        return a.getEstimate();
>    }
>    /**
>     * Get the value of the linear coefficient.
>     *
>     * @return the value of b for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getB() {
>        return b.getEstimate();
>    }
>    /**
>     * Get the value of the constant coefficient.
>     *
>     * @return the value of ac for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getC() {
>        return c.getEstimate();
>    }
>    /**
>     * Get the theoretical value of the model for some x.
>     * <p>The theoretical value is the value computed using
>     * the current state of the problem parameters.</p>
>     *
>     * Note the use of Hörner's method (synthetic division) for
> evaluating polynomials,
>     * (more efficient)
>     *
>     * @param x explanatory variable
>     * @return the theoretical value y = a x<sup>2</sup> + b x + c
>     */
>    public double theoreticalValue(double x) {
>        //System.out.println ("x = " + x + "  a.getEstimate() = " +
> a.getEstimate() + "  b.getEstimate() = " + b.getEstimate() + "
> c.getEstimate() = " + c.getEstimate());
>        return ( (a.getEstimate() * x + b.getEstimate() ) * x +
> c.getEstimate());
>    }
>    /**
>     * Get the partial derivative of the theoretical value
>     * of the model for some x.
>     * <p>The derivative is computed using
>     * the current state of the problem parameters.</p>
>     *
>     * @param x explanatory variable
>     * @param parameter estimated parameter (either a, b, or c)
>     * @return the partial derivative dy/dp
>     */
>    private double partial(double x, EstimatedParameter parameter) {
>        // since we know the only parameters are a, b and c in this
>        // class we simply use "==" for efficiency
>        if (parameter == a) {
>            return x * x;
>        } else if (parameter == b) {
>            return x;
>        } else {
>            return 1.0;
>        }
>    }
>    /** Internal measurements class.
>     * <p>The measurement is the y value for a fixed specified x.</p>
>     */
>    private class LocalMeasurement extends WeightedMeasurement {
>        static final long serialVersionUID = 1;
>        private final double x;
>        // constructor
>        public LocalMeasurement(double x, double y, double w) {
>            super(w, y);
>            this.x = x;
>        }
>        public double getTheoreticalValue() {
>            // the value is provided by the model for the local x
>            return theoreticalValue(x);
>        }
>        public double getPartial(EstimatedParameter parameter) {
>            // the value is provided by the model for the local x
>            return partial(x, parameter);
>        }
>    }
>  }

-- 
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.


Mime
View raw message