bugGNU Octave - Bugs: bug #50561, 350X slower code for eps() written...


bug #50561: 350X slower code for eps() written in C++ rather than an m-file

Submitted by:  Rik <rik5>
Submitted on:  Thu 16 Mar 2017 03:22:07 PM UTC  
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Performance
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  Closed
Release:  dev Operating System:  Any

Add a New Comment (Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


( Jump to the original submission)

Thu 16 Mar 2017 08:49:39 PM UTC, comment #8: 

All of the efficiency gains were due to removing the per-loop assignment of epsval to retval.  That was the change that was pushed.

After the cset was pushed, I tried a number of other optimizations but nothing made any recognizable change in execution times.

One can't return epsval because it is not of the polymorphic type octave_value which is the return type for all user-facing C++ functions.

Rik <rik5>
Project Administrator
Thu 16 Mar 2017 08:32:17 PM UTC, comment #7: 

What assignment above are you referring to?  Looking at the whole routine, I see now that

"retval = epsval"

is actually a large array assignment.  So clearly that was a huge reassignment of an array every time in the loop.

Did you try the things you mention before making the retval = epsval move?  After that retval=epsval move, maybe the more granular details you listed will now start to make a difference, improving things even further.

Also, if this retval is a large array, stack construct, does the C++ compiler optimize that out?  Or would it make sense to rather than assign retval = epsval, do things more direct with:

          return epsval;
              return epsval;
    return fill_matrix (args, std::numeric_limits<double>::epsilon (),


Dan Sebald <sebald>
Thu 16 Mar 2017 08:19:19 PM UTC, comment #6: 

I had already tried one optimization which was to declare min_double outside the loop.

double min_double = std::numeric_limits<double>::min ();

This didn't seem to make any difference.  Maybe because the value there is declared as a const and can be resolved at compile time.

I also tried substituting a function call

  epsval(i) = pow (2.0, -1074e0);

with a constant

  epsval(i) = 4.940656458412465e-324;

This made no difference either.  For test input I used

x = zeros (1e7, 1);

so the only path that was executed was the assignment above.  I was surprised by that one, but that's why you have to benchmark because human intuition can be faulty.

Rik <rik5>
Project Administrator
Thu 16 Mar 2017 06:27:51 PM UTC, comment #5: 

Let me add that even though a LUT looks to be a lot of memory, say, 1K or 2K, the machine code related to a pow() algorithm can be pretty big for an inefficient instruction base like x86.  Of course, in this application the pow() routines are present no matter what is done.  But if it were a case of substituting a 2K LUT for excluding a pow() routine in the app, it would be a fair trade.

Dan Sebald <sebald>
Thu 16 Mar 2017 06:11:20 PM UTC, comment #4: 

That "retval = epsval" was definitely redundant.  An fairly good optimizing compiler would take that out of the loop.

But there is more than that here.  In particular, this line:

  else if (val < std::numeric_limits<double>::min ())
      epsval(i) = pow (2.0, -1074e0);

I don't know how often that case gets called.  But pow(2.0,-1074e0) is a constant.  A really good optimizing compiler would take that outside the loop, but I doubt a fairly general compiler like linux gcc would do that.  (Only something like a C compiler for embedded chips would get that picky.)

Note that in your script, that is what is effectively done:

  retval(idx2) = pow2 (-1074);

The pow2 routine is called once, and then that value is assigned to all the appropriate elements.

So, in the C++ loop, the preference is to take that pow (2.0, -1074e0) prior to the loop by setting some stack variable and then re-using the stack variable.  Or even better, try to make it a static stack variable so that pow(2.0,-1074e0) resides in memory and then is simply a very fast machine instruction load.  The compiler might not like

static double smallpow = pow(2.0,-1074e0);

however, so you might need to be creative in how to generate that value at compile time.

Also, this pow2() and log2() function.  You may want to check where that is coming from.  It may be a noticeably faster implementation than the general routines:

pow (2.0, -1074e0);
octave::math::frexp (val, &expon);
std::pow (2.0,static_cast<double> (expon - 53));

For example, power of two computations can probably be done rather quick in integer math first with a single instruction bit shift.  (Above, there is a cast from integer to double, then two doubles for operators, then I'm sure general pow has much more work to do than simple bit shifts.)  If so, try substituting pow() with that routine specific to base-2 logarithmic computations.

Another possible little trick (who's performance probably depends on the nature of the data) would be to retain "expon" and "pow2(expon)" then reuse... hold on, why even do that?  Take a look at this particular instruction:

std::pow (2.0,static_cast<double> (expon - 53));

Because the left argument is fixed (2.0, i.e., 2), and I believe expon is an integer value with fairly limited range, -128 to 128, am I right?  The above line is ultimately a fairly short look-up table such that

std::pow (2.0,static_cast<double> (expon - 53));


pow2lut[expon - 53 + 128];

or whatever the proper offsets depending on the possible range of expon.

One could generate the LUT at startup, first time this routine is called, or during compilation by constructing the LUT

double pow2lut[256] = {

where the numbers could come from Octave, generated using full format precision.  The latter is a little dodgy for a platform portable program, probably only something I'd do for highly-controlled embedded code.  I think generating such a LUT the first time the routine is used though is fine.

Dan Sebald <sebald>
Thu 16 Mar 2017 04:03:27 PM UTC, comment #3: 

From earlier benchmarking, switching to using a pointer from data() will get about 10% improvement.  However, it is also less understandable code so I generally only made that change in liboctave.

However, I did check in an additional cset here which cleans up the eps implementation.  As long as we are using abs() which works for complex or real entries, there is no longer any need to call the C++ function fabs().  It didn't really save much in execution time, but it is more tidy.  See http://hg.savannah.gnu.org/hgweb/octave/rev/7a06a1a5a12b.

Rik <rik5>
Project Administrator
Thu 16 Mar 2017 03:45:25 PM UTC, comment #2: 

It would probably also be slightly faster to use double *pepsval = epsval.data(); and then pepsval[i] instead of using expsval(i) for indexing in the assignment inside the loop.  Or to use epsval.xelem(i).  Either of those avoids the check on the reference count.

John W. Eaton <jwe>
Project Administrator
Thu 16 Mar 2017 03:32:01 PM UTC, comment #1: 

I fixed the problem in this cset (http://hg.savannah.gnu.org/hgweb/octave/rev/aed81b364903).

The issue was that for the double case only, retval was being assigned the value of epsval for every iteration of the for loop.  The solution was just to bring it outside the loop.

Rik <rik5>
Project Administrator
Thu 16 Mar 2017 03:22:07 PM UTC, original submission:  

The eps is written in C++, and thus should be very fast.  However, it is easily beaten by an m-file implementation which is 350X faster and gets the same results.

What is going on here?  And does this explain why Octave seems to have been getting slower over the years?

The C++ code is in libinterp/corefcn/data.cc

          Array<double> x = arg0.array_value ();

          Array<double> epsval (x.dims ());

          for (octave_idx_type i = 0; i < x.numel (); i++)
              double val = ::fabs (x(i));
              if (octave::math::isnan (val) || octave::math::isinf (val))
                epsval(i) = lo_ieee_nan_value ();
              else if (val < std::numeric_limits<double>::min ())
                epsval(i) = pow (2.0, -1074e0);
                  int expon;
                  octave::math::frexp (val, &expon);
                  epsval(i) = std::pow (2.0,
                                        static_cast<double> (expon - 53));

              retval = epsval;

My m-file version

function [retval] = myeps (x)

  retval = abs (x);
  idx1 = isnan (retval) | isinf (retval);
  idx2 = retval < realmin;
  retval(idx1) = NaN;
  retval(idx2) = pow2 (-1074);
  idx3 = ! (idx1 | idx2);
  [~, exponent] = log2 (retval(idx3));
  retval(idx3) = pow2 (exponent - 53);


The results of testing

octave:45> x = randi (1e6, [1e5,1]);
octave:46> tic; y = eps (x); toc
Elapsed time is 3.175 seconds.
octave:47> tic; y2 = myeps (x); toc
Elapsed time is 0.00889206 seconds.
octave:48> 3.175/.00889206
ans =  357.06
octave:49> isequal (y,y2)
ans = 1
octave:50> diary off

Rik <rik5>
Project Administrator


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

Attach Files:

No files currently attached


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by sebald (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by rik5 (Submitted the item)

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.




    Follow 2 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-03-16 rik5 StatusNone => Fixed
        Open/ClosedOpen => Closed

    Back to the top

    Powered by Savane 3.5