bugGNU Scientific Library - Bugs: bug #40092, false position root finding...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #40092: false position root finding requires too many function evals

Submitter:  Patrick Alken <psa>
Submitted:  Mon 23 Sep 2013 03:49:24 PM UTC
   
 
Category:  Performance Severity:  3 - Normal
Operating System:  Status:  Postponed
Assigned to:  None Open/Closed:  Open
Release: 

Sat 05 Oct 2013 08:37:39 AM UTC, comment #5: 

The quick fix may be to check that x_linear differs
from x_left and x_right by some value. Some multiple of
GSL_DBL_EPSILON since required precision is not available
by design.

I do not think that additional tests of function values
may help since the derivative of the function can
be arbitrary large or small close to the root: 1e9*(x - 1)
or 1e-9*(x - 1).

> This clearly needs fixing but it requires design changes
> that will break backwards compatibility so we have
> unfortunately had to put this off to GSL version 2.x.


That is why I tried to raise discussion on gsl-help
instead of reporting a bug.

Max <nikulin>
Sat 05 Oct 2013 07:39:29 AM UTC, comment #4: 

Thank you for your further input.

As you say, the real problem is that the current design is not ideal because it does not make the desired precision available to the iteration routine.  And this can result in more iterations than are necessary.  

It appears that the residual test is intended to handle this situation but this also involves an extra function evaluation on each iteration because the solvers don't return function values.

This clearly needs fixing but it requires design changes that will break backwards compatibility so we have unfortunately had to put this off to GSL version 2.x.

Brian Gladman <gladman>
Group Member
Sat 05 Oct 2013 05:35:22 AM UTC, comment #3: 

I would like to stress the point from the original message.

The minimal step for root finding should not be much less
than the specified precision. The current design separates
the step and the test stages so the information
on the precision is not available when the new root
approximation is calculated. The result is more iteration
that for an efficient implementation of the same method.

Max <nikulin>
Sat 05 Oct 2013 05:20:09 AM UTC, comment #2: 

Brian,
Originally I faced the issue with a more cumbersome
function. I have tested the example with the polynomial
on linux (Intel Core 2 Duo, i686 kernel,
gcc (Ubuntu/Linaro 4.7.3-1ubuntu1) 4.7.3).

On linux the FPU works in 80 bit extended precision mode
by default. Windows set FPU to 64 bit double precision mode.
Due to precision loss, the test function evaluates to zero
at the point adjacent to the root that is representable
by the type double. It might be that the calculations
are performed by SSE2 instructions, not FPU.

I may reproduce the issue with my example by adding
-mpc64 option to gcc:
cc -mpc64 -g -Wall -W   gsl-falsepos.c  -lgsl -lgslcblas -lm -o gsl-falsepos

The bug is apparent if the test function is written
in the following way:
double fpol2 (double x, void *params)
{
   double *p = params;
   double rat = p[0]/p[1];
   return (x*p[2] + p[1])*(x + rat) - rat*p[2]*x;
}

Sorry, I forgot that 64 bit precision may lead to another
result.

(file #29298)

Max <nikulin>
Thu 03 Oct 2013 01:16:34 PM UTC, comment #1: 

Running the supplied program on Windows x64 gives:

# i N_f x_l-r x_i-r x_h-r x_h-x_l
  0     4 -2.000e-001  5.000e-012  5.000e-012  2.000e-001
  1     5  0.000e+000  0.000e+000  0.000e+000  0.000e+000
# Converged
# f(x_l) =  0.0000e+000, f(x_h) =  0.0000e+000

so I cannot reproduce the results set out in this report.

Brian Gladman <gladman>
Group Member
Mon 23 Sep 2013 03:49:24 PM UTC, original submission:  

From M.A.Nikulin@gmail.com:

Hi,

I have noticed that the GSL implementation of the false position method
for one-dimensional root finding may require twice more function
evaluations than the bisection method for an almost linear function that
should not be considered as a degenerate case.

Since the reference states:

> Its convergence is linear, but it is usually faster than bisection.


I think it worths a remark in the reference.

I suppose that it is a design flaw that user supplied tolerance is not
available in the iteration functions. The implementation of the Brent
method uses DBL_EPSILON, but in general user might request significantly
lower precision. As a result it deteriorates programs performance.

The problem appears if at a certain step the edge of a bracketing
interval almost coincides the root. Further linear interpolation gives
the same point (either due to rounding error or withing the required
precision). Function evaluation should be skipped for this point before
bisection step.

The following program demonstrates the issue. The bisection method
converges after ~50 function evaluations, the false position method
requires ~100 function evaluations for the same precision. I would
stress again that it is a well behaviored almost linear function.

Maxim Nikulin

/* Excessive function evaluations for an almost linear function
  * in the current implementation of the regula falsi method
  */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

double fpol2 (double x, void *params)
{
    double *p = params;
    return (x*p[2] + p[1])*x + p[0];
}

struct f_count_wrap_par {
    gsl_function *F;
    int n;
};

double f_count_wrap (double x, void *params)
{
    struct f_count_wrap_par *f = params;
    ++f->n;
    return GSL_FN_EVAL(f->F, x);
}

int main (void)
{
    int i, max_iter = 50;
    double x_i, x_l = 0.5, x_h = 1.1;
    double tol_abs = 5*GSL_DBL_EPSILON;
    double tol_rel = 5*GSL_DBL_EPSILON;
    gsl_root_fsolver *s;
    double par[] = { -1, 1, -1e-10 }; /* f(x) = -1e-10*x*x + (x - 1) */
    gsl_function F = {fpol2, par}, Fwrap;
    struct f_count_wrap_par wrap_par = {&F, 0};
    int status;
    double r = -2*par[0]/(sqrt(par[1]*par[1]-4*par[0]*par[2])+par[1]);

    Fwrap.function = f_count_wrap;
    Fwrap.params = &wrap_par;

    s = gsl_root_fsolver_alloc (gsl_root_fsolver_falsepos);
    /* s = gsl_root_fsolver_alloc (gsl_root_fsolver_bisection); */
    gsl_root_fsolver_set (s, &Fwrap, x_l, x_h);
    printf ("# i   N_f   x_l-r      x_i-r      x_h-r      x_h-x_l\n");
    for (i = 0, status = GSL_CONTINUE ;
          status == GSL_CONTINUE && i < max_iter ; ++i) {
       status = gsl_root_fsolver_iterate (s);
       x_i = gsl_root_fsolver_root (s);
       x_l = gsl_root_fsolver_x_lower (s);
       x_h = gsl_root_fsolver_x_upper (s);
       if (status != GSL_SUCCESS)
          break;
       status = gsl_root_test_interval (x_l, x_h, tol_abs, tol_rel);
       printf ("%3d %5d % .3e % .3e % .3e % .3e\n", i, wrap_par.n,
             x_l - r, x_i - r, x_h - r, x_h - x_l);
    }
    if (status == GSL_SUCCESS) {
       printf ("# Converged\n");
       printf ("# f(x_l) = % .4e, f(x_h) = % .4e\n",
             GSL_FN_EVAL(&F, x_l), GSL_FN_EVAL(&F, x_h));
    } else if (i == max_iter)
       printf ("# Max iteration number exceeded\n");
    gsl_root_fsolver_free (s);
    exit (EXIT_SUCCESS);
}

Output: iteration, number of function evaluations, shifts from the root
of the left interval edge, the current root estimation, and the right
edge, and inally the bracketing interval length

# i   N_f   x_l-r      x_i-r      x_h-r      x_h-x_l
   0     4 -2.000e-01  5.000e-12  5.000e-12  2.000e-01
   1     6 -1.000e-01  0.000e+00  0.000e+00  1.000e-01
   2     8 -5.000e-02  0.000e+00  0.000e+00  5.000e-02
   3    10 -2.500e-02  0.000e+00  0.000e+00  2.500e-02
   4    12 -1.250e-02  0.000e+00  0.000e+00  1.250e-02
   5    14 -6.250e-03  0.000e+00  0.000e+00  6.250e-03
   6    16 -3.125e-03  0.000e+00  0.000e+00  3.125e-03
   7    18 -1.563e-03  0.000e+00  0.000e+00  1.563e-03
   8    20 -7.813e-04  0.000e+00  0.000e+00  7.813e-04
   9    22 -3.906e-04  0.000e+00  0.000e+00  3.906e-04
  10    24 -1.953e-04  0.000e+00  0.000e+00  1.953e-04
  11    26 -9.766e-05  0.000e+00  0.000e+00  9.766e-05
  12    28 -4.883e-05  0.000e+00  0.000e+00  4.883e-05
  13    30 -2.441e-05  0.000e+00  0.000e+00  2.441e-05
  14    32 -1.221e-05  0.000e+00  0.000e+00  1.221e-05
  15    34 -6.104e-06  0.000e+00  0.000e+00  6.104e-06
  16    36 -3.052e-06  0.000e+00  0.000e+00  3.052e-06
  17    38 -1.526e-06  0.000e+00  0.000e+00  1.526e-06
  18    40 -7.629e-07  0.000e+00  0.000e+00  7.629e-07
  19    42 -3.815e-07  0.000e+00  0.000e+00  3.815e-07
  20    44 -1.907e-07  0.000e+00  0.000e+00  1.907e-07
  21    46 -9.537e-08  0.000e+00  0.000e+00  9.537e-08
  22    48 -4.768e-08  0.000e+00  0.000e+00  4.768e-08
  23    50 -2.384e-08  0.000e+00  0.000e+00  2.384e-08
  24    52 -1.192e-08  0.000e+00  0.000e+00  1.192e-08
  25    54 -5.960e-09  0.000e+00  0.000e+00  5.960e-09
  26    56 -2.980e-09  0.000e+00  0.000e+00  2.980e-09
  27    58 -1.490e-09  0.000e+00  0.000e+00  1.490e-09
  28    60 -7.451e-10  0.000e+00  0.000e+00  7.451e-10
  29    62 -3.725e-10  0.000e+00  0.000e+00  3.725e-10
  30    64 -1.863e-10  0.000e+00  0.000e+00  1.863e-10
  31    66 -9.313e-11  0.000e+00  0.000e+00  9.313e-11
  32    68 -4.657e-11  0.000e+00  0.000e+00  4.657e-11
  33    70 -2.328e-11  0.000e+00  0.000e+00  2.328e-11
  34    72 -1.164e-11  0.000e+00  0.000e+00  1.164e-11
  35    74 -5.821e-12  0.000e+00  0.000e+00  5.821e-12
  36    76 -2.910e-12  0.000e+00  0.000e+00  2.910e-12
  37    78 -1.455e-12  0.000e+00  0.000e+00  1.455e-12
  38    80 -7.276e-13  0.000e+00  0.000e+00  7.276e-13
  39    82 -3.637e-13  0.000e+00  0.000e+00  3.637e-13
  40    84 -1.819e-13  0.000e+00  0.000e+00  1.819e-13
  41    86 -9.104e-14  0.000e+00  0.000e+00  9.104e-14
  42    88 -4.552e-14  0.000e+00  0.000e+00  4.552e-14
  43    90 -2.265e-14  0.000e+00  0.000e+00  2.265e-14
  44    92 -1.132e-14  0.000e+00  0.000e+00  1.132e-14
  45    94 -5.773e-15  0.000e+00  0.000e+00  5.773e-15
  46    96 -2.887e-15  0.000e+00  0.000e+00  2.887e-15
  47    98 -1.332e-15  0.000e+00  0.000e+00  1.332e-15
# Converged
# f(x_l) = -1.3240e-15, f(x_h) =  8.2399e-18



Patrick Alken <psa>
Group administrator

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attached Files
file #29298:  gsl-falsepos64.c added by nikulin (2KiB - text/x-csrc - test program that works with -mpc64)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by nikulin (Updated the item)
  • -email is unavailable- added by gladman (Posted a comment)
  • -email is unavailable- added by psa (Submitted the item)
  •  

    Follow 3 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2013-10-05 gladman StatusWorks For Me Postponed
    2013-10-05 nikulin Attached File- Added gsl-falsepos64.c, #29298
    2013-10-03 gladman StatusNone Works For Me

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code