2
$\begingroup$

The problem occurs when setting the $\sigma$ of the GaussianFilter to zero and using the method Method -> "Gaussian".

The problem occurs here

GaussianFilter[RandomReal[1, {12, 101, 101}], {4 {.1, 0, 0}+ 1, {.1, 0, 0}},
  {2, 0, 0}, Padding -> "Periodic", Method -> "Gaussian"];

If we set the Method -> "Bessel" or if we set the $\sigma$ very low (something like 0.125) the computation is very fast.

Method Bessel fixes this

GaussianFilter[RandomReal[1, {12, 101, 101}], {4 {.1, 0, 0} + 1, {.1, 0, 0}},
  {2, 0, 0}, Padding -> "Periodic", Method -> "Bessel"];

Setting $\sigma$ to a non zero

GaussianFilter[RandomReal[1, {12, 101, 101}], {4 {.1, 0.125, 0.125}+ 1, {.1, 0.125, 0.125}},
  {2, 0, 0}, Padding -> "Periodic", Method -> "Gaussian"];

My Question is: Is this expected behavior or are we dealing with a bug? And can you explain why this is happening?

$\endgroup$
2
  • $\begingroup$ Please do not use the bugs tag until your observations have been confirmed by other users. $\endgroup$ Jan 31, 2017 at 10:27
  • $\begingroup$ Ok, was not aware of that. Thanks $\endgroup$ Jan 31, 2017 at 10:28

1 Answer 1

3
$\begingroup$

Timing

We can chart the effectiveness of the two methods with RepeatedTiming[]

range = {3, 5, 7, 10, 15, 20, 30, 40, 50};
bes = RepeatedTiming[
     GaussianFilter[
      RandomReal[
       1, {Round[#/3], #, #}], {4 {.1, 0, 0} + 1, {.1, 0, 0}}, {2, 0, 
       0}, Padding -> "Periodic", Method -> "Bessel"], 2][[1]] & /@ 
  range;
gauss = RepeatedTiming[
     GaussianFilter[
      RandomReal[
       1, {Round[#/3], #, #}], {4 {.1, 0, 0} + 1, {.1, 0, 0}}, {2, 0, 
       0}, Padding -> "Periodic", Method -> "Gaussian"], 2][[1]] & /@ 
  range;
ListLogPlot[{Transpose[{range, bes}], Transpose[{range, gauss}]}, 
 Joined -> True]

enter image description here

You can see that at i=50, the computation time climbs rapidly towards 10s and beyond. Compare this with a 2d version:

range2d = 10 range;
bes2d = RepeatedTiming[
     GaussianFilter[RandomReal[1, {#, #}], {{.1, 0, 0}}, {2, 0, 0}, 
      Padding -> "Periodic", Method -> "Bessel"], 2][[1]] & /@ range2d
gauss2d = RepeatedTiming[
     GaussianFilter[RandomReal[1, {#, #}], {{.1, 0, 0}}, {2, 0, 0}, 
      Padding -> "Periodic", Method -> "Gaussian"], 2][[1]] & /@ range2d
ListLogPlot[{Transpose[{range2d, bes2d}], Transpose[{range2d, gauss2d}]}, 
 Joined -> True]

enter image description here

If you increase the array sizes even further, the two lines become indistinguishable.

In other words, don't use Gaussian with (3+)d arrays.

So what's going on?

When we run without δ

GaussianFilter[
  RandomReal[1, {12, 101, 101}], {4 {.1, 0, 0} + 1, {.1, 0, 0}}, {2, 
   0, 0}, Padding -> "Periodic", Method -> "Gaussian"][[1, 1]]

Which gives a result like:

-0.558742 + (0.0931892 (0. + 14.7756 DiracDelta[0]^2))/DiracDelta[0]^2

This is an analytical result.

Whereas

GaussianFilter[
  RandomReal[
   1, {12, 101, 101}], {4 {.1, 10^-5, 10^-5} + 1, {.1, 10^-5, 
    10^-5}}, {2, 0, 0}, Padding -> "Periodic", 
  Method -> "Gaussian"][[1, 1]]

{-0.395175, 0.0122412,..., 0.469736, -0.102494}

Yields a numeric result. Something which produces underflow when the δ comes close to 0, meaning a numeric evaluation cannot be done when δ=0, which in higher dimensions causes a significant time increase.

This might mean that with a numeric evaluation you lose the added accuracy of the Gaussian method.

$\endgroup$
5
  • $\begingroup$ Thanks for you extensive research on the timing. I do want to add that the precision of the Gaussian is better $\endgroup$ Jan 31, 2017 at 12:22
  • $\begingroup$ @FrankMartin Then it's not that peculiar that it takes longer :). Your example takes my computer 6.72714s to run, which isn't excessive I think. $\endgroup$
    – Feyre
    Jan 31, 2017 at 12:37
  • $\begingroup$ It seems odd to me that when you set the $ \sigma $ to zero the computation time is longer compared when you set $\sigma$ to 0.125, I would expect the opposite. $\endgroup$ Jan 31, 2017 at 14:11
  • $\begingroup$ @FrankMartin And it is odd this just happens in 2-d. Nonetheless, maybe my edit to my answer gives some more insight. $\endgroup$
    – Feyre
    Jan 31, 2017 at 14:41
  • $\begingroup$ Yes that definitely helps, it gives some insight in how Mathematica is dealing with this internally. For now I'll use a small sigma as a work around. $\endgroup$ Jan 31, 2017 at 14:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.