I've retitled 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.

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 7e18 (those differences to Matlab were ~2e5 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.

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

Any updates on the 3.8.0 release?

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.

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

If the fix is indeed to rereverse 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 speedup the shape "same" which is currently implemented as computing the "full" shape and then clipping the borders.

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/octconvn.cc (around line 95, convolve_nd)
liboctave/cruft/blasxtra/dconv2.f (around line 73, daxpx call)

Any updates or roadmap on this one?

Yes, the higher dimension case makes thing more
difficult.
The way it is should be good enough.

The example keeps getting truncated.
I'll repost it since it came through successfully in the email.

I'm not sure that it is possible to account for the "order" problem the way Ndimensional 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:

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.

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:

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:

With hg id > "97ed9dd479ab tip @", I'm seeing 7 conv2 test failures
libinterp/corefcn/conv2.cctst ......................... PASS 40/47 FAIL 7
I've attached my fntests.log
(file #29429)

Already ahead of you :) those are already fixed
http://hg.savannah.gnu.org/hgweb/octave/rev/7903ae3b657d

Also, did you notice that in some of the
failed tests the array dimensions do not match?
Could be a problem with the test?

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:

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.

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/octconvn.cc) and exactly inside the following loops:
The 'b + mb*jb' argument should be 'b + mb*(nbjb1)' in fact convolution is always a sum of products of this type: a(u1, u2, u3) * b(x1u1, x2u2, x3u3).
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)

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.

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.

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.

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.

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.29912184192067e005 1.28215490056687e005 1.29375784844437e005
1.20384206842329e005 1.60724929372802e005 1.55126096841249e005
1.10716816072338e005 1.79941888939234e005 1.73947117905957e005
1.67516463603148e005 2.38806994258701e005 2.04159342070756e005
2.49735609200875e005 2.98445469683259e005 2.23170170254184e005
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.
