[jira] Created: (MATH-352) Jacobian rank determination in LevenbergMarquardtOptimizer is not numerically robust

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

[jira] Created: (MATH-352) Jacobian rank determination in LevenbergMarquardtOptimizer is not numerically robust

JIRA jira@apache.org
Jacobian rank determination in LevenbergMarquardtOptimizer is not numerically robust
------------------------------------------------------------------------------------

                 Key: MATH-352
                 URL: https://issues.apache.org/jira/browse/MATH-352
             Project: Commons Math
          Issue Type: Bug
    Affects Versions: 2.0
         Environment: commons-math 2.0
            Reporter: Gene Gorokhovsky
            Priority: Critical
             Fix For: 2.1


LevenbergMarquardtOptimizer is designed to handle singular jacobians,  i.e. situations when some of the fitted parameters depend on each other. The check for that condition is in LevenbergMarquardtOptimizer.qrDecomposition uses precise comparison to 0.

    if (ak2 == 0 ) {
                rank = k;
                return;
        }

A correct check would be comparison with a small epsilon. Hard coded 2.2204e-16 is used elsewhere in the same file for similar purpose.

final double QR_RANK_EPS = Math.ulp(1d); //2.220446049250313E-16
....
    if (ak2  < QR_RANK_EPS) {
                rank = k;
                return;
        }

Current exact equality check is not tolerant of the real world poorly conditioned situations. For example I am trying to fit a cylinder into sample 3d points. Although theoretically cylinder has only 5 independent variables, derivatives for optimizing function (signed distance) for such minimal parametrization are complicated and it  it much easier to work with a 7 variable parametrization (3 for axis direction, 3 for axis origin and 1 for radius). This naturally results in rank-deficient jacobian, but because of the numeric errors the actual ak2 values for the dependent rows ( I am seeing values of 1e-18 and less), rank handling code does not kick in.
Keeping these tiny values around then leads to huge corrections for the corresponding very slowly changing parameters, and consequently to numeric errors and instabilities. I have noticed the problem because tiny shift in the initial guess (on the order of 1e-12 in the axis component and origins) resulted in significantly different finally converged answers (origins and radii differing by as much as 0.02) which I tracked to loss of precision due to numeric error with root cause described above.
Providing a cutoff as suggested fixes the issue. After the fix, small perturbations in the initial guess had practically no effect to the converged result - as expected from a robust algorithm.


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

Reply | Threaded
Open this post in threaded view
|

[jira] Updated: (MATH-352) Jacobian rank determination in LevenbergMarquardtOptimizer is not numerically robust

JIRA jira@apache.org

     [ https://issues.apache.org/jira/browse/MATH-352?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]

Phil Steitz updated MATH-352:
-----------------------------

    Fix Version/s:     (was: 2.1)
                   2.2

> Jacobian rank determination in LevenbergMarquardtOptimizer is not numerically robust
> ------------------------------------------------------------------------------------
>
>                 Key: MATH-352
>                 URL: https://issues.apache.org/jira/browse/MATH-352
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 2.0
>         Environment: commons-math 2.0
>            Reporter: Gene Gorokhovsky
>            Priority: Critical
>             Fix For: 2.2
>
>
> LevenbergMarquardtOptimizer is designed to handle singular jacobians,  i.e. situations when some of the fitted parameters depend on each other. The check for that condition is in LevenbergMarquardtOptimizer.qrDecomposition uses precise comparison to 0.
>     if (ak2 == 0 ) {
>                 rank = k;
>                 return;
>         }
> A correct check would be comparison with a small epsilon. Hard coded 2.2204e-16 is used elsewhere in the same file for similar purpose.
> final double QR_RANK_EPS = Math.ulp(1d); //2.220446049250313E-16
> ....
>     if (ak2  < QR_RANK_EPS) {
>                 rank = k;
>                 return;
>         }
> Current exact equality check is not tolerant of the real world poorly conditioned situations. For example I am trying to fit a cylinder into sample 3d points. Although theoretically cylinder has only 5 independent variables, derivatives for optimizing function (signed distance) for such minimal parametrization are complicated and it  it much easier to work with a 7 variable parametrization (3 for axis direction, 3 for axis origin and 1 for radius). This naturally results in rank-deficient jacobian, but because of the numeric errors the actual ak2 values for the dependent rows ( I am seeing values of 1e-18 and less), rank handling code does not kick in.
> Keeping these tiny values around then leads to huge corrections for the corresponding very slowly changing parameters, and consequently to numeric errors and instabilities. I have noticed the problem because tiny shift in the initial guess (on the order of 1e-12 in the axis component and origins) resulted in significantly different finally converged answers (origins and radii differing by as much as 0.02) which I tracked to loss of precision due to numeric error with root cause described above.
> Providing a cutoff as suggested fixes the issue. After the fix, small perturbations in the initial guess had practically no effect to the converged result - as expected from a robust algorithm.

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

Reply | Threaded
Open this post in threaded view
|

[jira] Resolved: (MATH-352) Jacobian rank determination in LevenbergMarquardtOptimizer is not numerically robust

JIRA jira@apache.org
In reply to this post by JIRA jira@apache.org

     [ https://issues.apache.org/jira/browse/MATH-352?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]

Luc Maisonobe resolved MATH-352.
--------------------------------

    Resolution: Fixed

Fixed in subversion repository as of r951864.
A setQRRankingThreshold has been added as proposed in [http://markmail.org/message/p2j76cnwsyehl7u6].
Thanks for reporting the issue.

> Jacobian rank determination in LevenbergMarquardtOptimizer is not numerically robust
> ------------------------------------------------------------------------------------
>
>                 Key: MATH-352
>                 URL: https://issues.apache.org/jira/browse/MATH-352
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 2.0
>         Environment: commons-math 2.0
>            Reporter: Gene Gorokhovsky
>            Priority: Critical
>             Fix For: 2.2
>
>
> LevenbergMarquardtOptimizer is designed to handle singular jacobians,  i.e. situations when some of the fitted parameters depend on each other. The check for that condition is in LevenbergMarquardtOptimizer.qrDecomposition uses precise comparison to 0.
>     if (ak2 == 0 ) {
>                 rank = k;
>                 return;
>         }
> A correct check would be comparison with a small epsilon. Hard coded 2.2204e-16 is used elsewhere in the same file for similar purpose.
> final double QR_RANK_EPS = Math.ulp(1d); //2.220446049250313E-16
> ....
>     if (ak2  < QR_RANK_EPS) {
>                 rank = k;
>                 return;
>         }
> Current exact equality check is not tolerant of the real world poorly conditioned situations. For example I am trying to fit a cylinder into sample 3d points. Although theoretically cylinder has only 5 independent variables, derivatives for optimizing function (signed distance) for such minimal parametrization are complicated and it  it much easier to work with a 7 variable parametrization (3 for axis direction, 3 for axis origin and 1 for radius). This naturally results in rank-deficient jacobian, but because of the numeric errors the actual ak2 values for the dependent rows ( I am seeing values of 1e-18 and less), rank handling code does not kick in.
> Keeping these tiny values around then leads to huge corrections for the corresponding very slowly changing parameters, and consequently to numeric errors and instabilities. I have noticed the problem because tiny shift in the initial guess (on the order of 1e-12 in the axis component and origins) resulted in significantly different finally converged answers (origins and radii differing by as much as 0.02) which I tracked to loss of precision due to numeric error with root cause described above.
> Providing a cutoff as suggested fixes the issue. After the fix, small perturbations in the initial guess had practically no effect to the converged result - as expected from a robust algorithm.

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