How can I collect the sum of any number of gradients without using ReplaceRepeated[]?

Question

Recently, I have run into the problem where I need to collect the the sum of any number of gradients. I have defined a function grad[] to represent the gradient, and a typical question is like:

how to collect

 grad[f]-grad[g]+grad[h]-grad[p]

into

grad[f-g+h-p]

Gradient can be collected this way since it is a linear operator. And so do other linear operators. I know I can write a binary rule and use ReplaceRepeated[] like:

grad[f]-grad[g]+grad[h]-grad[p]//.grad[a_]+grad[b_]:>grad[a+b]

But the problem is, I am not satisfied with the algorithm efficiency of ReplaceRepeated[]. What I want is a patterned rule which does the trick with using ReplaceAll[] once rather than repeatedly.

Any idea how to achieve this?

Answer 1

linearCombine from Pushing Mathematica's FullSimplify to a global complexity minimum does the job:

ClearAll[linearCombine];
Module[{f}, SetAttributes[f, NumericFunction];
  constantQ = NumericQ[# /.
      s_Symbol /; MemberQ[Attributes[s], Constant] :> f[0]] &];
linearCombine[e_Plus, head_, constants_List : {}] := 
  Block[constants, SetAttributes[#, Constant] & /@ constants;
   head@Replace[e, {c_?constantQ*head[a_] :> c*a, head[a_] :> a}, 1] /;
     MatchQ[e, _[((_?constantQ)*_head | _head) ..]]];

linearCombine[grad[f] - grad[g] + grad[h] - grad[p], grad]

(*  grad[f - g + h - p]  *)

Answer 2

You might try TagSet, something like this:

grad /: Times[-1, grad[x_]] := grad[-x]
grad /: Plus[ps___, grad[x_], qs___, grad[y_], rs___] := Plus[grad[x + y], ps, qs, rs]

The definition of Times might need to be extended depending on what sorts of expressions you expect to encounter. Also, some care might need to be taken to avoid nesting of grad. The definition for Plus might be overkill--I didn't do any real testing.