Let's say I'm trying to write my own gradient operator function.

If I know the function is going to be $\mathbb{R}^2\rightarrow\mathbb{R}$, I could write

grad[f_Function] := { D[f[x,y],x] , D[f[x,y],y] }

But what if I want to make grad work for functions $\mathbb{R}^n\rightarrow\mathbb{R}$?

I could do something like this:

grad[f_Function] := Module[{xargs = Table[Symbol["x"<>ToString[i]],{i,1,n}]},
    Table[D[f[xargs],xargs[[j]]], {j, 1, n}]

But then this only works with f[{x1,x2,x3,x4}] instead of f[x1,x2,x3,x4]. Moreover, how should grad know what n should be?


4 Answers 4


Update: Alternatively, you can define an Operator

SetAttributes[gradOp, HoldAll]
gradOp[f_[x___]] := Through@Through[(Derivative @@@ IdentityMatrix[Length@{x}])@f][x]
gradOp@f[1, 2, 1, t] // TeXForm

$ \left\{f^{(1,0,0,0)}(1,2,1,t),f^{(0,1,0,0)}(1,2,1,t),f^{(0,0,1,0)}(1,2,1,t),f^{(0,0,0,1)}(1,2,1,t)\right\} $


gradOp2 = Module[{vars = Array[\[FormalX], Length@#]}, 
    Through[Grad[#[[0]] @@ vars, vars][[;; , 0]][## & @@ #]]] &;

gradOp2@f[1, 2, s, t] // TeXForm

$ \left\{f^{(1,0,0,0)}(1,2,s,t), f^{(0,1,0,0)}(1,2,s,t),f^{(0,0,1,0)}(1,2,s,t), f^{(0,0,0,1)}(1,2,s,t)\right\} $

Original answer:

gradF = Grad[#, List @@ #] &;

gradF@f[ w, x, y, z]

{(f^(1,0,0,0))[w,x,y,z], (f^(0,1,0,0))[w,x,y,z], (f^(0,0,1,0))[w,x,y,z], (f^(0,0,0,1))[w,x,y,z]}


$\left\{f^{(1,0,0,0)}(w,x,y,z),f^{(0,1,0,0)}(w,x,y,z),f^{(0,0,1,0)} (w,x,y,z),f^{(0,0,0,1)}(w,x,y,z)\right\} $

Note: In case you need it in other contexts, arg need be a List, for Length[arg] to work. For example,

Length[qwertyiop[u, v,  w, x, y, z]]



Here's one way:

z = Array[x, 3];
D[f[z] /. List -> Sequence, #] & /@ z

Define your variables as x[1], x[2], ... (as many as you wish, and call them z). The first line uses the Array function to do this succinctly. Then apply the D function to f[z]. The List->Sequence solves the brackets issue (i.e., replaces f[{x1,x2,x3,x4}] with f[x1,x2,x3,x4]) and the use of map (the funny symbols /@) means it will be applied no matter how many variables there are.


You could try this:

(* This method only works with functions with an explicit list of arguments,
     not with functions that use slots (#) *)
f = Function[{x, y, z}, x y z]; 

grad[funct_Function] := With[{
    n = Length[funct[[1]]] (* find the number of arguments*)
        (* construct the derivatives *)
        Derivative @@@ 
          Reverse @ Sort @ Permutations[Prepend[ConstantArray[0, n - 1], 1]]

Evaluate the gradient:


Apply the gradient to arguments:

Through[grad[f][1, 1, 2]]

If you want this version of grad to work with slots, you could try to look into the function and find the largest instance of Slot[i_] by using Cases. Obviously that doesn't work with SlotSequence (##):

f = Function[#1 + #3]
Max[Cases[f, Slot[i_Integer] :> i, \[Infinity]]]

For explicit function arguments, you can write it like this:

grad[f_[args__]] := Table[D[f[args], a], {a, {args}}]

This does exactly what you want. For example:

grad[f[x, y]]
(* {(f^(1,0))[x,y],(f^(0,1))[x,y]} *)


grad[f[x, y, z, w]]
(* {(f^(1,0,0,0))[x,y,z,w],(f^(0,1,0,0))[x,y,z,w],(f^(0,0,1,0))[x,y,z,w],(f^(0,0,0,1))[x,y,z,w]} *)

The trick is the args__ pattern, which means that args becomes the Sequence of arguments, no matter how long it is (but there must be at least one argument). You can convert this Sequence into a List by writing {args}. This you can then use, e.g., in a Table the way I did. For completeness: the "dimension" of the argument function is Length[{args}].


If you want a gradient operator instead, you could write

gradOp[f_[args__]] := With[{dim = Length[{args}]},
  Table[Derivative[Sequence@@UnitVector[dim, n]][f][args], {n, dim}]

such that

gradOp@f[12, y]
(* {(f^(1,0))[12,y],(f^(0,1))[12,y]} *)

You can even make this operator work with pure functions, as long as they don't contain a SlotSequence. First, make sure the above definition does not evaluate for pure function arguments, using /;.

gradOp[f_[args__]] := With[{dim = Length[{args}]},
   Table[Derivative[Sequence @@ UnitVector[dim, n]][f][args], {n, dim}]
] /; FreeQ[f[args], Slot[_] | SlotSequence[___]]

Then, define the case of a pure-function argument.

gradOp[f_Function] := Module[{C},
  (* extract all slot names in f *)
  With[{slots = Cases[f, Slot[i_] :> i, \[Infinity]]},
   (* replace the slots by variables *)
   With[{expr = f /. Thread[Thread[Slot[slots]] -> Thread[C[slots]]] /. Function -> Identity},
    (* construct a new pure function... *)
      (* ...which represents the gradient of f *)
      Table[D[expr, a], {a, Thread[C[slots]]}] /. Thread[Thread[C[slots]] -> Thread[Slot[slots]]]
] /; FreeQ[f, SlotSequence[___]]

With this gradOp works on both, explicit and pure functions.

(* {f'[x]} *)

gradOp[#1 - #2^2 &]
(* {1, -2 #2} & *)

gradOp[#1 - #2^2 &][x, y]
(* {1, -2 y} *)

It even keeps the slot names:

gradOp[Log[#xyz] + #test^2 &]
(* {1/#xyz, 2 #test}& *)

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.