I have this very simple function where I'm trying to find the Jacobian, but all my different solutions do not agree with the computed numerical difference Jacobian.
The function I'm looking at is:
$$\mathbf{v} = 0.5 \operatorname{sinc}\left(\frac{|\mathbf{x}|}{2}\right) \: \mathbf{x} $$
And the jacobian is, according to my calculations:
$$\frac{\partial \mathbf{v}}{\partial \mathbf{x}} = 0.5 \left(\mathbf{x}\mathbf{x}^T \frac{2}{|\mathbf{x}|^3}(\frac{|\mathbf{x}|}{2}\cos(\frac{|\mathbf{x}|}{2}) - \sin(\frac{|\mathbf{x}|}{2}) + \operatorname{sinc}(\frac{|\mathbf{x}|}{2})\:\mathbf{I}\right) $$
but somehow it's wrong. A short Python script to test:
import numpy as np
def f(x):
return 0.5 * x * np.sinc(0.5 * np.linalg.norm(x))
x0 = np.array([0.0288708, 0.0189388, 0.0823266])
J_numDiff = np.empty((3, 3), dtype=float)
# compute numerical difference jacobian with perturbation d
d = 1e-8
J_numDiff[:, 0] = (f(x0 + np.array([d, 0, 0])) - f(x0))/d
J_numDiff[:, 1] = (f(x0 + np.array([0, d, 0])) - f(x0))/d
J_numDiff[:, 2] = (f(x0 + np.array([0, 0, d])) - f(x0))/d
# compute analytical jacobian with helpers
norm = np.linalg.norm(x0)
halfnorm = 0.5 * norm
x0 = x0[:, np.newaxis]
xxT = x0 * x0.T
J_analytic = 0.5 * (xxT * 2 / (norm**3) * (np.cos(halfnorm) * halfnorm - np.sin(halfnorm))
+ np.eye(3) * np.sinc(halfnorm))
# difference is not zero
diff = J_numDiff - J_analytic
print(diff)
Numpy is using $$\textrm{sinc}x=\frac{\sin\pi x}{\pi x}$$ So create a function
and call it instead of
np.sinc. You need to do this both in the definition of $f$ as well when you calculate the analytic formula. The differences will be on the order of 1e-10, which are numerical errors