Applying Operators inside complex expressions to argument with Through

Question

I am developing some functions to help me work with Quantum Mechanics Momentum Operators. Here, for simplicity, I represent them with these simple functions because the actions of the operators themselves is not the problem.

f1[x_] := c1*x
f2[x_] := c2*x
f3[x_] := c3*x
f4[x_] := c4*x

These operators/functions appear alone or as compositions in algebraic expressions and inside vectors. Something like this:

vec = {f1@*f3 - f2@*f4, f3 + 2 f4 + 4, const1 + Exp[x] - f1};

Now I need to apply this vector to an argument [a], and I'm pretty sure that the main function to achieve this goal is Through. The correct end result should be:

{a c1 c3 - a c2 c4, a c3 + 2 a c4 + 4, const1 + Exp[x] - a c1}

Firs step:

vec2 = Through[vec[a]]
{(f1@*f3 - f2@*f4)[a], (4 + f3 + 2 f4)[a], (const1 + E^x - f1)[a]}

This first step works as expected, feeding this argument [a] to each component of the vector.

Now I need this to keep happening inside the expressions of each component of the vector, so that the operators f1,f2,.. get evaluated with the argument [a]. Constants and anything that is not an operator should remain unchanged throughout the process. To do this I Map Through to that output, but now I get mixed results:

Map[Through, vec2]
{a c1 c3 + (-(f2@*f4))[a], a c3 + 4[a] + (2 f4)[a], 
 const1[a] + (E^x)[a] + (-f1)[a]}

It works fine when the argument hits the operator directly, with operators alone or with composition, but fails otherwise:

1. When anything goes with the operators/compositions: (-(f2@*f4))[a] or (2 f4)[a] or (-f1)[a], because the operator doesn't get evaluated.

2. It feeds the argument to everything, not just operators: 4[a], const1[a] + (E^x)[a].

I have tried playing with the second argument of Through (Through[expr,h] performs the transformation wherever h occurs in the head of expr.), with other functions different to Map and with different levelspecs for Map and other functions, but honestly I don't think sharing those results would help.

I will add that, if I apply Through again to a relevant part like:

Through[(-(f2@*f4))[a]]
a c2 c4 (-1)[a]

The operator composition does get evaluated but it stills try to also evaluate (-1)[a]. Another thing is that I am not sure how to do this other than by hand. This makes me think that a proper recursive/multilevel mapping with the correct head argument of Through (or other kind of discrimination) would do the trick. And what drives me crazy is that I am almost sure that it will be a single line of code! Any help would be greatly appreciated. Thanks a lot!

Edit: Explicit list of operators & ReplaceAll approach

The more I study it, the more robust and sound seem to me the implementation proposed by @Lukas Lang based on the use of a specific head to identify operators. What I love of this approach is that it now looks easy to add special behaviours to handle more and more operations, like Dot products, Cross products of operators or combined actions of operators on |j,m> or TensorProduct[|j1,m1>,|j2,m2>] spaces, etc. In the following days I will try to develop it further, at least for problems involving Angular Momentum and Addition of Angular Momentum.

I'm pretty convinced that it's not as sound as @Lujas Lang approach, but for educational purposes I would like to share another possibility that I was considering. It involves maintaining an explicit list of operators and then use a function that use a simple pattern matching set of rules to "apply" operators to arguments when needed:

ClearAll[f1, f2, f3, f4, opApply, opList]
opList = f3 | f2 | f1 | f4;
f1[x_] := c1*x
f2[x_] := c2*x
f3[x_] := c3*x
f4[x_] := c4*x
opApply[list_List[args__]] := opApply /@ Through[list[args]]
opApply[head_[args__]] := 
 ReplaceAll[Expand[head], {op1_@*op2_ /; (MatchQ[op1, opList] && MatchQ[op2, opList]) :> op1[op2[args]], op1_ /; MatchQ[op1, opList] :> op1[args]}]

vec = {f1@*f3 - f2@*f4, f3 + 2 f4 + 4, const1 + Exp[x] - f1};
vec[a]
% // opApply

Output:

{f1@*f3 - f2@*f4, 4 + f3 + 2 f4, const1 + E^x - f1}[a]
{a c1 c3 - a c2 c4, 4 + a c3 + 2 a c4, -a c1 + const1 + E^x}

It's not well tested, it's just to show the idea.

On the one hand, it's cumbersome in it's own way: you have to maintain the list of operators and the list of rules may get a little bit crazy or difficult to read if you want to handle many operations or situations.

On the other hand, it looks that have the potential to let you work freely with operators in a very natural notation and just use opApply when you need to as a special form of Through if you wish.

Since I am learning Wolfram Language, I would love to hear some opinions on this approach. Advantages? Disadvantages? Could it work good in complex situations? Should it be avoided for any reason?

Answer 1

I think the easiest way to do this is to simply implement all rules of the operator arithmetic explicitly, rather than trying to be too clever about it. For example, the + operator in vec is not in fact the standard $+:\mathbb{C}\times\mathbb{C}\rightarrow\mathbb{C}$, but rather $+:\mathbb{(C\rightarrow C)}\times\mathbb{(C\rightarrow C)}\rightarrow\mathbb{(C\rightarrow C)}$. So in some sense, the rules anyway need to be redefined on a case-by-case basis.

Here's how I would do it: First, choose a head operator that designates that an expression should follow our operator arithmetic. We define what happens if this operator is applied to a list of arguments by implementing special cases for the common operations (using the tricks from the question):

operator[op : (List | Plus | Times)[___]][args___] := Through[(operator /@ op)[args]]
operator[op : (Composition | RightComposition)[___]][args___] := (operator /@ op)[args]

We then define a fallback rule for constants:

operator[const_][_] := const

Finally, we define our custom operators:

operator[f1][x_] := c1 x
operator[f2][x_] := c2 x
operator[f3][x_] := c3 x
operator[f4][x_] := c4 x

Now, the composite operator vec:

vec = operator[{f1@*f3 - f2@*f4, f3 + 2 f4 + 4, const1 + Exp[x] - f1}]
(* operator[{f1@*f3 - f2@*f4, 4 + f3 + 2 f4, const1 + E^x - f1}] *)

Applying it to a:

vec[a]
(* {a c1 c3 - a c2 c4, 4 + a c3 + 2 a c4, -a c1 + const1 + E^x} *)

You'll note that what we're doing here is essentially the multilevel/recursive mapping you were suggesting in your question, just with a helper head for readability and ease of use. Also, you'll probably need to add some more rules for arithmetic operations, e.g. Power.

Further remarks

One potential pitfall that comes to mind is your distinction of multiplication and composition: Usually in QM, operators are always linear, and so a*b is used to unambiguously refer to the composition of operators a and b (since (a*b)[c] must be a[b[c]] and not a[c]*b[c], as the latter would not be linear). If you choose to adopt that notation, constants would instead be operators multiplying their argument with themselves, so e.g.: 2[c] is 2*c. This notation has the advantage that a^3 and by extension Exp[a] can be directly interpreted as a*a*a and 1+a+a*a/2+…, respectively. Finally, (1+a)[c] is then c+a[c] which is again linear, rather than 1+a[c], which is not.

An implementation of this would look something like this:

Clear@operator

operator[op : (List | Plus)[___]][args___] := Through[(operator /@ op)[args]]
operator[op : (Composition | Times | NonCommutativeMultiply)[___]][args___] := (operator /@ Composition @@ op)[args]
operator[const_][x_] := const x

operator[f1][x_] := c1 x
operator[f2][x_] := c2 x
operator[f3][x_] := c3 x
operator[f4][x_] := c4 x

vec = operator[{f1 ** f3 - f2 ** f4, f3 + 2 f4 + 4, const1 + Exp[x] - f1}]
(* operator[{f1 ** f3 - f2 ** f4, 4 + f3 + 2 f4, const1 + E^x - f1}] *)

vec[a]
(* {a c1 c3 - a c2 c4, 4 a + a c3 + 2 a c4, -a c1 + a const1 + a E^x} *)

The end result is of course slightly different in this case, but as noted, this way, the vec operator is linear, whereas in your interpretation it is not (which might or might not be desired of course). Also note the use of ** (NonCommuativeMultiply) instead of *, since of course operator multiplication/composition is in general not commutative. Nevertheless, I implemented the second rule in such a way that it supports *, **, and @* for increased flexibility, while interpreting all of them as @*.

Remarks on edit

Here are some remarks on the edit to the question, and the alternate approach proposed there:

• First of all, the second rule can be made quite a bit more readable by replacing op1_ /; MatchQ[op1, opList] by the more direct op1: opList (for the documentation, see Pattern):

  opApply[head_[args__]] := 
   ReplaceAll[Expand[head], {(op1:opList)@*(op2:opList) :> op1[op2[args]], op1:opList :> op1[args]}]
  

• In general, I prefer the approach shown in this answer, because it is (at least in my experience) often more readable and leverages Mathematica's strength in pattern matching better by separating the different rules into separate downvalues. However, sometimes an approach based on ReplaceAll is also the better solution, especially in situations involving tight evaluation control. And of course, it's also a matter of personal preference.

• Furthermore, I think the approach is not so different in it's core, with the two main differences being:

  • The use of ReplaceAll instead of individual rules for all expressions: This enables easier traversal of deeper expressions, since we don't need to implement rules for all levels. Of course, this gives us less control, but it's also less work if we don't need the control. We could emulate this for the operator approach as follows (untested):

    operator[h_[parts___]][args__] := (operator[h][args]) @@ (operator[#][args] & /@ {parts})
    operator[atom_][__] := atom
    

    The first rule essentially makes sure operator is "applied" to every part of an expression, while the second rule effectively terminates the traversal of the expression once the bottom is reached. All other rules would then correspond to the actual replacement rules

  • The use of a function that is applied to a full expression as opposed to a helper head that is wrapped around the operator part. We can also emulate this for the operator approach:

    opApply[op_[args__]]:=operator[op][args]
    

    Clearly this is cheating in some sense, but it does somehow combine the ease of use of the opApply approach with the ease of definition of the operator approach.