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

Submitter:  Alberto Andreotti <albertoandreotti>
Submitted:  Fri 21 Jun 2013 09:05:08 PM UTC
   
 
Category:  Octave Function Severity:  2 - Minor
Priority:  3 - Low Item Group:  Inaccurate Result
Status:  Confirmed Assigned to:  None
Originator Name:  AlbertoAndreotti Open/Closed:  * Open
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 09 Aug 2023 08:42:55 PM UTC, comment #26: 

updating that this behavior (errors ~ 10-20eps in convn) is still present in 8.2.0

Nicholas Jankowski <nrjank>
Group Member
Mon 21 Nov 2016 11:35:13 PM UTC, 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>
Group administrator
Mon 21 Nov 2016 10:05:34 PM UTC, 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 02 Dec 2013 02:16:58 PM UTC, 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>
Group Member
Mon 02 Dec 2013 10:28:02 AM UTC, comment #22: 

Any updates on the 3.8.0 release?

Anonymous
Tue 12 Nov 2013 10:20:19 PM UTC, 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>
Group administrator
Tue 12 Nov 2013 09:24:18 PM UTC, comment #20: 

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

Anonymous
Tue 12 Nov 2013 09:00:02 PM UTC, 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


changeset:   10946:1094868ca10d
user:        Jaroslav Hajek <highegg@gmail.com>
date:        Tue Sep 07 12:23:01 2010 +0200
summary:     fix bugs in inner convolution


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


changeset:   14596:ad1d7a53b1ca
branch:      stable
parent:      14590:6250e1232c9c
user:        Sunghyun Cho <sodomau@gmail.com>
date:        Tue Apr 17 12:47:34 2012 +0100
summary:     perform convolution properly independently of matrix orientation (bug #34893)


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>
Group Member
Tue 12 Nov 2013 07:46:12 PM UTC, 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>
Group administrator
Tue 12 Nov 2013 07:35:45 PM UTC, comment #17: 

Any updates or roadmap on this one?

Anonymous
Sun 03 Nov 2013 06:06:52 PM UTC, comment #16: 

Yes, the higher dimension case makes thing more
difficult.

The way it is should be good enough.

Michael Godfrey <godfrey>
Group Member
Sun 03 Nov 2013 05:37:11 PM UTC, comment #15: 

The example keeps getting truncated.

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


octave> a = rand (1, 27);
octave> b = rand (1, 27);
octave> c1 = sum (a .* b);
octave> c2 = sum (a(end:-1:1) .* b(end:-1:1));
octave> c1 - c2
ans =  1.7764e-15


Rik <rik5>
Group administrator
Sun 03 Nov 2013 05:15:42 PM UTC, 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:


octave> a = rand (1, 27);
octave> b = rand (1, 27);
octave> c1 = sum (a .* b);
octave> c2 = sum (a(end:-1:1) .* b(end:-1:1));
octave> c1 - c2
ans =  1.7764e-15
-verbatim

Carnë Draug <carandraug>
Group Member
Sun 03 Nov 2013 04:44:42 PM UTC, 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>
Group Member
Mon 21 Oct 2013 01:22:24 PM UTC, 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:


octave> a = rand (1, 27);
octave> b = rand (1, 27);
octave> c1 = sum (a .* b);
octave> c2 = sum (a(end:-1:1) .* b(end:-1:1));
octave> c1 - c2
ans =  1.7764e-15
-verbatim

Carnë Draug <carandraug>
Group Member
Sun 20 Oct 2013 07:29:42 PM UTC, 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:


a = rand (10, 10, 10);
b = rand (3, 3, 3);
full  = convn (a, b, "full");
valid = convn (a, b, "valid");
d = full(3:10,3:10,3:10) - valid; # these should be exactly the same
max (d(:))
min (d(:))


I get before the patch:


octave> max (d(:))
ans =  1.4640
octave> min (d(:))
ans = -1.4087


and after the patch:


octave> max(d(:))
ans =    3.5527e-15
octave> min(d(:))
ans =   -2.6645e-15


Carnë Draug <carandraug>
Group Member
Sun 20 Oct 2013 07:04:06 PM UTC, 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>
Group Member
Sun 20 Oct 2013 04:11:20 PM UTC, comment #9: 

Already ahead of you :) those are already fixed

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

Carnë Draug <carandraug>
Group Member
Sun 20 Oct 2013 03:59:29 PM UTC, 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>
Group Member
Sun 20 Oct 2013 03:55:09 PM UTC, 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:


a = rand (10, 10, 10);
b = rand (3, 3, 3);
full  = convn (a, b, "full");
valid = convn (a, b, "valid");
d = full(3:10,3:10,3:10) - valid; # these should be exactly the same
max (d(:))
min (d(:))


I get before the patch:


octave> max (d(:))
ans =  1.4640
octave> min (d(:))
ans = -1.4087

and after the patch:

+verbatim+
octave> max(d(:))
ans =    3.5527e-15
octave> min(d(:))
ans =   -2.6645e-15


Carnë Draug <carandraug>
Group Member
Sat 19 Oct 2013 06:37:39 AM UTC, 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>
Group Member
Mon 12 Aug 2013 10:25:55 PM UTC, 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:


for (octave_idx_type ja = 0; ja < na - nb + 1; ja++)
  for (octave_idx_type jb = 0; jb < nb; jb++)
    convolve_nd<T, R> (a + ma*(ja + jb), ad, acd, b + mb*jb, bd, bcd,
                       c + ldc*ja, ccd, nd-1, inner);


The 'b + mb*jb' 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:


octave:1> v = [1,2]
v =

   1   2

octave:2> vvv = reshape(v,1,1,2)
vvv =

ans(:,:,1) =  1
ans(:,:,2) =  2

octave:3> convn(vvv, vvv, 'valid')
ans =  5
octave:4> convn(vvv, vvv)
ans =

ans(:,:,1) =  1
ans(:,:,2) =  4
ans(:,:,3) =  4


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 19 Jul 2013 03:02:21 PM UTC, 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>
Group administrator
Fri 19 Jul 2013 12:26:27 PM UTC, 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 19 Jul 2013 06:14:07 AM UTC, comment #2: 

Convolution has been no end of trouble...

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


conv2full = convn (zzz, hhh, 'full');
conv2same = convn (zzz, hhh, '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>
Group administrator
Mon 24 Jun 2013 07:28:41 PM UTC, 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):


y1 = convn (zzz, hhh, "full")(8:12, 8:12, 5);
y2 = convn (zzz, hhh, "same")(4:8, 4:8, 3);
y3 = convn (zzz, hhh, "valid");
assert (y2, y1);
assert (y3, y1);


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>
Group Member
Fri 21 Jun 2013 09:05:08 PM UTC, 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 Files:
   
   
Comment:
   

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

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

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

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 15 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2016-11-21 rik5 Severity3 - Normal 2 - Minor
        Priority5 - Normal 3 - Low
        Summaryconvn(x,y,&quot;valid&quot;) produces results that are different from Matlab convn over arrays of greater than 2 dimensions accumulates errors on the order of 20*eps
    2013-12-14 carandraug Dependencies- bugs #40895 is dependent
    2013-10-20 bpabbott Attached File- Added fntests.log, #29429
    2013-10-20 carandraug Open/ClosedClosed Open
    2013-10-20 carandraug StatusFixed Confirmed
    2013-10-19 carandraug StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2013-08-12 mcecchetti Attached File- Added convn.diff, #28829
    2013-06-24 mtmiller Item 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
    2013-06-21 albertoandreotti Attached File- Added octave_convn.m, #28380

    Back to the top

    Powered by Savane 3.13-cf05.
    Corresponding source code