Incorrect second derivative of sinc nearby zero

[petvana]

As I have already mentioned in julia#37227, the function sinc is numerically unstable nearby zero. The sinc function is especially important for this value because its smoothness is the reason to be utilized, but the second derivative is not computed correctly.

using ForwardDiff
f(x) = ForwardDiff.derivative(y -> ForwardDiff.derivative(sinc, y), x)
for i in [1:5:89..., 90:20:330...] 
  x = 0.1^i
  println("f(", x, ") = ", f(x))
end

Results:

f(0.1) = -3.193029835964822
f(1.0000000000000004e-6) = -3.28955078125
f(1.0000000000000006e-11) = -2.097152e6
f(1.0000000000000008e-16) = 1.8014398509481984e16
f(1.0000000000000012e-21) = 0.0
f(1.0000000000000015e-26) = -2.658455991569832e36
f(1.0000000000000016e-31) = 3.4253944624943037e46
f(1.000000000000002e-36) = 1.9615942923083377e56
f(1.0000000000000023e-41) = 0.0
f(1.0000000000000026e-46) = 0.0
f(1.0000000000000028e-51) = 3.7299242730734e86
f(1.0000000000000031e-56) = 0.0
f(1.0000000000000033e-61) = 0.0
f(1.0000000000000037e-66) = -6.304320991423117e116
f(1.000000000000004e-71) = 0.0
f(1.0000000000000042e-76) = -3.488825876618913e136
f(1.0000000000000045e-81) = -2.8418048210060247e160
f(1.0000000000000047e-86) = Inf
f(1.000000000000005e-90) = Inf
f(1.0000000000000061e-110) = Inf
f(1.0000000000000073e-130) = Inf
f(1.0000000000000083e-150) = Inf
f(1.0000000000000094e-170) = NaN
f(1.0000000000000106e-190) = NaN
f(1.0000000000000117e-210) = NaN
f(1.0000000000000127e-230) = NaN
f(1.0000000000000139e-250) = NaN
f(1.000000000000015e-270) = NaN
f(1.0000000000000162e-290) = NaN
f(1.0e-310) = NaN
f(0.0) = 0.0

The second derivative is well defined and the limit at zero is

-(pi ^ 2) / 3 = -3.289868133696453

Is it possible to provide a specific rule for this case? It would be nice to be able to use it directly for NLP in JuMP.

[andreasnoack]

I'm not sure how to easily fix this. It might require a good implementation of the derivative.

[stevengj]

The next version of Julia should have an accurate sinc derivative implementation (cosc) (JuliaLang/julia#37273) that is supported in ChainRules (JuliaDiff/ChainRules.jl#255). If you AD the new cosc implementation it should be accurate near zero, I think, since the new code is just implemented as a polynomial near zero with no catastrophic cancellations.

[andreasnoack]

@petvana Then we just need a definition for sinc in https://github.com/JuliaDiff/DiffRules.jl/blob/master/src/rules.jl and then we can benefit from the improved cosc once it has made it into a release.

[stevengj]

I think this issue can be closed now. With the original code above, I now get:

f(0.1) = -3.1930298359648535
f(1.0000000000000004e-6) = -3.2898681336867113
f(1.0000000000000006e-11) = -3.2898681336964524
f(1.0000000000000008e-16) = -3.2898681336964524
f(1.0000000000000012e-21) = -3.2898681336964524
f(1.0000000000000015e-26) = -3.2898681336964524
f(1.0000000000000016e-31) = -3.2898681336964524
f(1.000000000000002e-36) = -3.2898681336964524
f(1.0000000000000023e-41) = -3.2898681336964524
f(1.0000000000000026e-46) = -3.2898681336964524
f(1.0000000000000028e-51) = -3.2898681336964524
f(1.0000000000000031e-56) = -3.2898681336964524
f(1.0000000000000033e-61) = -3.2898681336964524
f(1.0000000000000037e-66) = -3.2898681336964524
f(1.000000000000004e-71) = -3.2898681336964524
f(1.0000000000000042e-76) = -3.2898681336964524
f(1.0000000000000045e-81) = -3.2898681336964524
f(1.0000000000000047e-86) = -3.2898681336964524
f(1.000000000000005e-90) = -3.2898681336964524
f(1.0000000000000061e-110) = -3.2898681336964524
f(1.0000000000000073e-130) = -3.2898681336964524
f(1.0000000000000083e-150) = -3.2898681336964524
f(1.0000000000000094e-170) = -3.2898681336964524
f(1.0000000000000106e-190) = -3.2898681336964524
f(1.0000000000000117e-210) = -3.2898681336964524
f(1.0000000000000127e-230) = -3.2898681336964524
f(1.0000000000000139e-250) = -3.2898681336964524
f(1.000000000000015e-270) = -3.2898681336964524
f(1.0000000000000162e-290) = -3.2898681336964524
f(1.0e-310) = -3.2898681336964524
f(0.0) = -3.2898681336964524