# Speed up derivative evaluation

I'm trying to calculate the normal vectors and the tangent vectors at the discrete points of a suface. For example:

f[x_, y_] := (Cos[x] Sin[y])/2. (* The surface.*)
fx[a_, b_] :=
Module[{x, y}, D[f[x, y], x] /. {x -> a, y -> b}]
fy[a_, b_] :=
Module[{x, y}, D[f[x, y], y] /. {x -> a, y -> b}]
normalVector[a_, b_] :=
Module[{x, y},
{-D[f[x, y], x], -D[f[x, y], y], 1} /. {x -> a, y -> b}]

tangentVector[x_, y_, θ_] :=
{Cos[θ], Sin[θ],
fx[x, y] Cos[θ] + fy[x, y] Sin[θ]};

xr = yr = 1;
n = 20;
hx = 2 xr/n;
hy = 2 yr/n;
hθ = 2 Pi/n;
AbsoluteTiming[
Table[
nV = normalVector[x, y];
tV = tangentVector[x, y, θ];
{nV, tV},
{x, -xr, xr - hx, hx}, {y, -yr, yr - hy, hy},
{θ, -Pi, Pi - hθ, hθ}
];
]


gives (vectors shown only at one point)

The problem is that the calculation is slow. The desired discrete point should be 200 per axis. I think the time is mostly consumed by derivative evaluation. So how could I speed up the calculation ? One choice is to transform the code into C++. However the surface function might come from interpolation. I want to stay in Mathematica.

• You know that the gradient can be computed in one go, right? D[f[x, y], {{x, y}}] – J. M.'s ennui May 26 '15 at 4:04
• Nice tip for speedup. – novice May 26 '15 at 4:09

As simple as it sound, use Set instead of delayed assignment. i.e., replace your := in the function definitions with =.

First, make sure everything is safe by clearing x,y,a,b as:

ClearAll[x,y,a,b]


Then, define your functions as:

f[x_, y_] = (Cos[x] Sin[y])/2.; (*The surface.*)

fx[a_, b_] = Module[{x, y}, D[f[x, y], x] /. {x -> a, y -> b}];
fy[a_, b_] = Module[{x, y}, D[f[x, y], y] /. {x -> a, y -> b}];
normalVector[a_, b_] =
Module[{x,
y}, {-D[f[x, y], x], -D[f[x, y], y], 1} /. {x -> a, y -> b}];


I got 5x improvement in speed just by doing that.

This works because the function f[x,y] is fully defined and all the derivatives can be obtained symbolically beforehand. What is happening with the delayed assignment, is basically having D[f[x,y],x] being calculated each time a call is made for fx[a,b] is made. Repetitive evaluation get cashed, but apparently still not good enough in this case.

Also, I am not sure about your C++ comment, but a great idea would be to compile your code from Mathematica. Define f[x,y] as above and then use:

f[x_, y_] = (Cos[x] Sin[y])/2.; (*The surface.*)
fx = Compile[{{a, _Real}, {b, _Real}},
Module[{x, y}, D[f[x, y], x] /. {x -> a, y -> b}]];
fy = Compile[{{a, _Real}, {b, _Real}},
Module[{x, y}, D[f[x, y], y] /. {x -> a, y -> b}]];


Update: After experimenting with Compile, there is no performance gain obtained. Only the use of Set (=) instead of the delayed assignment (:=).

• @Bichoy, great improvement ! another 15x speedup. Does compile support interpolation functions? – novice May 26 '15 at 3:07
• I am not sure :) You have to try ... You should also read the documentation for Compile, you can set the compilation target to "C", set the RuntimeAttributes to {Listable} and enable parallelization for even much more improvement... Just check Compile documentation. – Bichoy May 26 '15 at 3:11
• @novice Sorry, I had a syntax error in the compiled version (when copying) please recheck that. – Bichoy May 26 '15 at 3:14
• @@Bichoy, I think Compile fails for functions containing D. The compiled version is slower. – novice May 26 '15 at 4:59
• @novice, @bichoy, always worth checking this list of compilable functions before looking into Compile – dr.blochwave May 26 '15 at 6:45

Here's one speed-up: Use approximate machine real inputs. Outer is a bit faster than Table here.

Clear[x, y];
f[x_, y_] := (Cos[x] Sin[y])/2 (*The surface.*)
df[x_, y_] = D[f[x, y], {{x, y}}];

normalVector[a_, b_] := Join[-df[a, b], {1}];
tangentVector[x_, y_, θ_] := Join[#, {df[x, y].#}] &@{Cos[θ], Sin[θ]};

xr = yr = 1;
n = 20.;        (* < this makes the Range's below machine reals *)
hx = 2 xr/n;
hy = 2 yr/n;
hθ = 2 Pi/n;
RepeatedTiming[
foo = Outer[{normalVector[#1, #2], tangentVector[##]} &,
Range[-xr, xr - hx, hx],
Range[-yr, yr - hy, hy],
Range[-Pi, Pi - hθ, hθ]];]
(*   {0.11, Null}   *)


If you value speed over readability, here's this, which takes some advantage of the vectorization of Mathematica functions. More savings may be possible.

RepeatedTiming[
glurg = Outer[
Transpose@{
ConstantArray[normalVector[#1, #2], Round[2 Pi/hθ]],
Transpose @ tangentVector[#1, #2, Range[-Pi, Pi - hθ, hθ]]} &,
Range[-xr, xr - hx, hx], Range[-yr, yr - hy, hy]];]
(*   {0.019, Null}   *)

foo == glurg
(*  True  *)

• @@Michael E2, the idea of vectorization is great. I get 20x speedup. I am curious about where does the time-saving come from? What if normalVector is also a function of θ? – novice May 27 '15 at 1:14
• @novice, for many built-in math functions, f[array] is handled in the math library (in low-level code) and is several times faster than, say, Map[f, array, {-1}]. See mathematica.stackexchange.com/a/21863 for more tips. As for normalVector, it would depend on how it is a function of θ. – Michael E2 May 27 '15 at 1:25

Just for fun:

f[x_, y_] := Cos[x] Sin[y]/2
grd[u_, v_] := Grad[z - f[x, y], {x, y, z}] /. {x -> u, y -> v};
tg[x0_, y0_, {a_, b_}] :=
With[{ru = {1, 0, D[f[x, y], x] /. {x -> x0, y -> y0}},
rv = {0, 1, D[f[x, y], y] /. {x -> x0, y -> y0}}}, a ru + b rv]
tn[x_, y_, {s_, t_}] :=
With[{pt = {x, y, f[x, y]}},
Graphics3D[{Point[pt], PointSize[0.04], Arrowheads[0.02], Red,
Arrow[{pt, pt + Normalize[grd[x, y]]}], Blue,
Arrow[{pt, pt + Normalize[tg[x, y, {s, t}]]}]}]]
tab = Flatten[
Table[Show[
Plot3D[f[x, y], {x, 0, 2 Pi}, {y, 0, 2 Pi},
PlotRange -> {-1, 1.5}, PlotStyle -> Opacity[0.5],
ColorFunction -> "DarkRainbow"],
tn[#, #, {Cos[t], Sin[t]}] &@pnt], {pnt, 0, 2 Pi, 0.3}, {t, 0,
2 Pi, 0.3}], 1];


An individual plot does not take long. The following gif has been cut down from tab` to display: