Possible problem with custom derivative of jnp.sinc. at x=0

[jecampagne]

Hello,
I may do not anderstand correctly the custom_jvp:

jax/jax/_src/numpy/ufuncs.py

Lines 2940 to 2953 in 15024ba

	@partial(custom_jvp, nondiff_argnums=(0,))
	def _sinc_maclaurin(k, x):
	# compute the kth derivative of x -> sin(x)/x evaluated at zero (since we
	# compute the monomial term in the jvp rule)
	# TODO(mattjj): see https://github.com/jax-ml/jax/issues/10750
	if k % 2:
	return x * 0
	else:
	return x * 0 + _lax_const(x, (-1) ** (k // 2) / (k + 1))
	@_sinc_maclaurin.defjvp
	def _sinc_maclaurin_jvp(k, primals, tangents):
	(x,), (t,) = primals, tangents
	return _sinc_maclaurin(k, x), _sinc_maclaurin(k + 1, x) * t

But looking at the successive derivative at x=0

sincp = grad(jnp.sinc)
sincp2 = grad(sincp)
sincp3 = grad(sincp2)
sincp4 = grad(sincp3)
sincp(0.),sincp2(0.),sincp3(0.),sincp4(0.)

gives

(Array(0., dtype=float64, weak_type=True),
 Array(-3.28986813, dtype=float64, weak_type=True),
 Array(0., dtype=float64, weak_type=True),
 Array(19.48181821, dtype=float64, weak_type=True))

But clearly in Mathematica

sinc(Pi x)^{(p)}(x=0) = If[Mod[p, 2] == 1, 0, Pi^p (-1)^(p/2)/(p + 1)!]

which will gives for the sincp2(0.), sincp4(0.).

-1.64493, 0.811742

And if one instead uses

 If[Mod[p, 2] == 1, 0, Pi^p (-1)^(p/2)/(p + 1)]}

then one recovers the current JAX sinc behavior as

-3.28987, 19.4818

So I wander if someones forget to use the factorial "(p+1)!" (ie. Gamma[p+2]) and use simply ("p+1").

If I am true this is a bug in the current version. and I propse

@partial(custom_jvp, nondiff_argnums=(0,)) 
 def _sinc_maclaurin(k, x): 
   # compute the kth derivative of x -> sin(x)/x evaluated at zero (since we 
   # compute the monomial term in the jvp rule) 
   # TODO(mattjj): see https://github.com/jax-ml/jax/issues/10750 
   if k % 2: 
     return x * 0 
   else: 
     return x * 0 + _lax_const(x, (-1) ** (k // 2) / jax.numpy.exp(jax.lax.lgamma(k+2.))

If not then there is some I do not undersrtand.

[jakevdp]

I don't understand the mathematica result you quote: it just seems wrong.

For example, here is a table of sinc derivatives at x = 0 that agrees with the formula and output values used in JAX, and disagrees with the formula and values you mention from mathematica: https://calculus.subwiki.org/wiki/Sinc_function#Higher_derivatives_at_zero Note that this table is for $\sin(x)/x$, while JAX uses the convention $\sin(\pi x)/(\pi x)$ so you need to multiply the values in the table by a factor of $\pi^k$:

>>> np.pi ** 2 * (-1 / 3)
-3.289868133696453

>>> np.pi ** 4 * (1 / 5)
19.481818206800483

Checking another source, wolfram alpha returns values that match this table, and thus agrees with JAX: https://www.wolframalpha.com/input?i=table+d%5En%2Fdx%5En+sinc%28x%29+for+n+%3D+1+...+4

[jecampagne]

Well in Mathematuca

Series[Sinc[x  Pi], {x, 0, 6}]

and

Table[{p, If[Mod[p, 2] == 1, 0, Pi^p (-1)^(p/2)/(p + 1)!]}, {p, 
  Range[5]}]

While

Table[{p, If[Mod[p, 2] == 1, 0, Pi^p (-1)^(p/2)/(p + 1)]}, {p, 
   Range[5]}] // N

gives

{{1., 0.}, {2., -3.28987}, {3., 0.}, {4., 19.4818}, {5., 0.}}

[jakevdp]

I would double-check the step where you go from the series expansion to the first table. When you start taking nth derivatives of the series with respect to x and evaluate at x = 0, the factorial terms cancel out.

[jecampagne]

Ho! my bad :( Yes of course you're right. I was focusing in something else. Sorry.

[jakevdp]

No worries, thanks for flagging the potential issue!