bugGNU Octave - Bugs: bug #39314, convn over arrays of greater than...

 
 

bug #39314: convn over arrays of greater than 2 dimensions accumulates errors on the order of 20*eps

Submitted by:  Alberto Andreotti <albertoandreotti>
Submitted on:  Fri Jun 21 21:05:08 2013  
 
Category: Octave FunctionSeverity: 2 - Minor
Priority: 3 - LowItem Group: Inaccurate Result
Status: ConfirmedAssigned to: None
Originator Name: AlbertoAndreottiOpen/Closed: Open
Release: devOperating System: GNU/Linux

Add a New Comment (Rich MarkupRich Markup):
   

You are not logged in

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

 

(Jump to the original submission Jump to the original submission)

Mon Nov 21 23:35:13 2016, comment #25:

I've re-titled this report to reflect the remaining issue. I've also lowered the severity and priority. If anyone wants to work on this they can, but errors on the order of tens of eps are not surprising nor unworkable.

Rik <rik5>
Project Administrator
Mon Nov 21 22:05:34 2016, comment #24:

The original issue of this bug report has been mostly (but not yet fully) solved in Octave 4.2.0.

My results on the code in comment #0 are now all below 7e-18 (those differences to Matlab were ~2e-5 originally).

But the slight imprecision in convolutions (of dimension > 2) is still present. With the code in comment #11 I still get values around 10 to 20 eps.

Hartmut <hardy>
Mon Dec 2 14:16:58 2013, comment #23:

Yes, there are updates on the 3.8 release. Please refer to the maintainer's list and the Octave wiki.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Mon Dec 2 10:28:02 2013, comment #22:

Any updates on the 3.8.0 release?

Anonymous
Tue Nov 12 22:20:19 2013, comment #21:

I agree. I think instead of reversing we just need to supply the correct offset.

The fix will be part of 3.8.0 which is due for release in a few weeks.

Rik <rik5>
Project Administrator
Tue Nov 12 21:24:18 2013, comment #20:

Ok, sounds good. When will the fix be released? As part of what version?

Anonymous
Tue Nov 12 21:00:02 2013, comment #19:

If the fix is indeed to re-reverse the order of the summation, then the fix may involve reverting the following changesets:

http://hg.savannah.gnu.org/hgweb/octave/rev/1094868ca10d

http://hg.savannah.gnu.org/hgweb/octave/rev/ad1d7a53b1ca

I'm unsure why reverting it was required. I may be recalling wrong, but I think simply adding an extra offset when passing the pointers to convolve_nd should be enough. This would avoid having a separate routine for the valid shape, and possibly speed-up the shape "same" which is currently implemented as computing the "full" shape and then clipping the borders.

Carnë Draug <carandraug>
Project Member
Tue Nov 12 19:46:12 2013, comment #18:

Octave has been patched so it produces correct results to within a few tens of eps. This is far better than producing incorrect values which it was before. It should be possible to eliminate even the small difference of a few tens of eps, but this is rather less urgent. The problem has, I believe, been correctly identified by Carne in that the order of summation differs when doing inner versus outer convolution.

For those interested, the relevant files are

liboctave/numeric/oct-convn.cc (around line 95, convolve_nd)
liboctave/cruft/blas-xtra/dconv2.f (around line 73, daxpx call)

Rik <rik5>
Project Administrator
Tue Nov 12 19:35:45 2013, comment #17:

Any updates or roadmap on this one?

Anonymous
Sun Nov 3 18:06:52 2013, comment #16:

Yes, the higher dimension case makes thing more
difficult.

The way it is should be good enough.

Michael Godfrey <godfrey>
Sun Nov 3 17:37:11 2013, comment #15:

The example keeps getting truncated.

I'll re-post it since it came through successfully in the e-mail.

Rik <rik5>
Project Administrator
Sun Nov 3 17:15:42 2013, comment #14:

I'm not sure that it is possible to account for the "order" problem the way N-dimensional convolution is being done. convn reduces the problem to a series of 2 dimensional convolutions, whose results are added.

You will also notice that that there's no precision errors when limited to 2 dimensions.

I see now that I have missed the markup on my previous comment. The example was:

Carnë Draug <carandraug>
Project Member
Sun Nov 3 16:44:42 2013, comment #13:

Has comment #12 been verified? Why would the order
of summation on average change the accuracy? High
accuracy summation algorithms take account of the
"order" problem by, in effect, computing partial sums.
But, I am sure this is not done, or appropriate, here.

In any case, it should be checked that the required
tolerance is within the expected rounding error.

Michael Godfrey <godfrey>
Mon Oct 21 13:22:24 2013, comment #12:

I think the cause of the different values comes from the order the values are added. Looking at the code for inner convolution and comparing to other, it seems that it's being done in reverse order so machine precision errors get in. Kind of like in the following example:

Carnë Draug <carandraug>
Project Member
Sun Oct 20 19:29:42 2013, comment #11:

When I said on comment #9 "Already ahead of you :) those are already fixed" I meant the failing tests due to wrong array dimensions mentioned on comment #8. cset 7903ae3b657d "convn: fix incorrect tests" only fixed the tests that were wrong.

But the bug about precision as mentioned on comment #7 still holds and that's why the tests are failing (this is not a regression, the tests are failing because the tests are new but the bug has always been there).

I see now that I messed up the markup on comment #7 so half of the the comment is gone. Here is it again.


Compare the following before and after the patch:

I get before the patch:

and after the patch:

Carnë Draug <carandraug>
Project Member
Sun Oct 20 19:04:06 2013, comment #10:

With hg id -> "97ed9dd479ab tip @", I'm seeing 7 conv2 test failures

libinterp/corefcn/conv2.cc-tst ......................... PASS 40/47 FAIL 7

I've attached my fntests.log

(file #29429)

Ben Abbott <bpabbott>
Project Member
Sun Oct 20 16:11:20 2013, comment #9:

Already ahead of you :) those are already fixed

http://hg.savannah.gnu.org/hgweb/octave/rev/7903ae3b657d

Carnë Draug <carandraug>
Project Member
Sun Oct 20 15:59:29 2013, comment #8:

Also, did you notice that in some of the
failed tests the array dimensions do not match?

Could be a problem with the test?

Michael Godfrey <godfrey>
Sun Oct 20 15:55:09 2013, comment #7:

I'm reopening this bug. I thought the tests were passing but apparently I messes up somewhere. With the change suggested we can't see the bug with integers anymore. Comparing the results to before the patch, the fix is correct but somewhere along the way some precision is lost (but close to machine precision). Compare the following before and after the patch:

I get before the patch:

Carnë Draug <carandraug>
Project Member
Sat Oct 19 06:37:39 2013, comment #6:

Thank you for your patch. I have pushed the change you suggested and added tests to convn. It seems to pass all of them.

From what I understood, this patch will pass a pointer to the last of each dimension of B. I noticed that cset 1094868ca10d inverted the order of accessing elements for valid shape only, and later ad1d7a53b1ca tried to fix that. I don't know Fortran to understand it well enough, just reporting this here in case someone thinks it may make sense to change the fix to somewhere else.

Carnë Draug <carandraug>
Project Member
Mon Aug 12 22:25:55 2013, comment #5:

Indeed the problem can be reproduced with integer values too and is present only for arrays with more than 2 dimensions. The bug is inside the convolve_nd function (liboctave/numeric/oct-convn.cc) and exactly inside the following loops:

The 'b + mbjb' argument should be 'b + mb(nb-jb-1)' in fact convolution is always a sum of products of this type: a(u1, u2, u3) * b(x1-u1, x2-u2, x3-u3).

The following octave example is straightforward:

The correct result for convn(vvv, vvv, 'valid') should be vvvv(1,1,2) that is 4 instead we get the same result as we computed dot(v, v).

I attached a patch with the correct implementation. I tested it with the zzz and hhh as input and it works fine.

(file #28829)

Marco Cecchetti <mcecchetti>
Fri Jul 19 15:02:21 2013, comment #4:

There's no need to upload the Octave versions of the matrices; I can generate those on my side. I think it would be useful to have Matlab's version because I want an independent idea of what the calculation should yield. Whoever originally coded the convolution routine thought they had it correct and I don't want to "think" I'm doing the right thing and get it wrong as well.

Rik <rik5>
Project Administrator
Fri Jul 19 12:26:27 2013, comment #3:

Hi Rik, guys,

I think you can compute them from the data that I already uploaded. But if you want me to put both matlab/octave versions, let me know.
Now, I'm at work, but tomorrow I can do it if necessary.

Cheers,

Alberto.

Alberto Andreotti <albertoandreotti>
Fri Jul 19 06:14:07 2013, comment #2:

Convolution has been no end of trouble...

For reference, can the original reporter also upload the matrices for "full" and "same"?

I took a look at the code and we are doing something different for the 'valid' parameter. It may be that we are not calling the Fortran correctly, or maybe the underlying BLAS call isn't right. At least superficially 'valid' seems to work, but the test cases I have are with integer values for the two matrices A and B.

Rik <rik5>
Project Administrator
Mon Jun 24 19:28:41 2013, comment #1:

Thank you for your bug report. Confirmed, also present in the development version of Octave. Note that the error you are seeing is present only in the "valid" mode of the convn function (and possibly conv2).

The error can also be observed this way (using your provided dataset for this example):

This shows that it is an error in how the "valid" convolution subset is computed. The "full" and "same" operations give the same result (within a much tighter tolerance) as the conv2d matrix you have supplied that was generated by Matlab.

Mike Miller <mtmiller>
Project Administrator
Fri Jun 21 21:05:08 2013, original submission:

Hi,

the convn() function is returning results that are different from the ones returned on Matlab. The error is small, but causes problems in some applications.
To reproduce it, do,

>load octave_convn.m
>convn(zzz, hhh, 'valid') - conv2d


zzz, hhh, are the input matrices, and conv2d the output of convn(zzz,hhh, 'valid') in Matlab. All this stuff comes in the file, and was generated in Matlab.
This, when run on Matlab, will return a matrix of zeros, and something like this(error)when run in Octave,

1.29912184192067e-005 1.28215490056687e-005 1.29375784844437e-005
1.20384206842329e-005 1.60724929372802e-005 1.55126096841249e-005
1.10716816072338e-005 1.79941888939234e-005 1.73947117905957e-005
1.67516463603148e-005 2.38806994258701e-005 2.04159342070756e-005
2.49735609200875e-005 2.98445469683259e-005 2.23170170254184e-005

the above response is of course truncated(two columns missing).
This represents the difference between Octave's computation and Matlabs.
Although the error is small, it is enough to prevent the convn() function to be used in many applications. This was tested in Ubuntu and in Windows 7 with the same results.

Alberto Andreotti <albertoandreotti>

 

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

Attach File(s):
   
   
Comment:
   

Attached Files
file #29429:  fntests.log added by bpabbott (313kB - application/octet-stream)
file #28829:  convn.diff added by mcecchetti (1kB - text/x-patch)
file #28380:  octave_convn.m added by albertoandreotti (8kB - application/octet-stream - the file was generated in Matlab.)

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -unavailable- added by hardy (Posted a comment)
  • -unavailable- added by jordigh (Posted a comment)
  • -unavailable- added by bpabbott (Updated the item)
  • -unavailable- added by godfrey (Posted a comment)
  • -unavailable- added by mcecchetti (Updated the item)
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by mtmiller (Posted a comment)
  • -unavailable- added by albertoandreotti (Submitted the item)
  •  

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

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 15 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Mon Nov 21 23:35:13 2016rik5Severity3 - Normal=>2 - Minor
      Priority5 - Normal=>3 - Low
      Summaryconvn(x,y,"valid") produces results that are different from Matlab=>convn over arrays of greater than 2 dimensions accumulates errors on the order of 20*eps
    Sat Dec 14 17:56:04 2013carandraugDependencies-=>bugs #40895 is dependent
    Sun Oct 20 19:04:06 2013bpabbottAttached File-=>Added fntests.log, #29429
    Sun Oct 20 15:55:36 2013carandraugOpen/ClosedClosed=>Open
    Sun Oct 20 15:55:09 2013carandraugStatusFixed=>Confirmed
    Sat Oct 19 06:37:39 2013carandraugStatusConfirmed=>Fixed
      Open/ClosedOpen=>Closed
    Mon Aug 12 22:25:55 2013mcecchettiAttached File-=>Added convn.diff, #28829
    Mon Jun 24 19:28:41 2013mtmillerItem GroupIncorrect Result=>Inaccurate Result
      StatusNone=>Confirmed
      Release3.6.4=>dev
      SummaryThe convn() function produces results that are different from Matlab.=>convn(x,y,"valid") produces results that are different from Matlab
    Fri Jun 21 21:05:08 2013albertoandreottiAttached File-=>Added octave_convn.m, #28380

    Back to the top


    Powered by Savane 3.1-cleanup