# Define multilinear operator

In Maple it is possible to define a multilinear, orderless operator with:

define(A,'orderless',multilinear):


I can't seem to find the equivalent in Mathematica.

I would like it to do something like:

A[x,x]/.x->0+b x+c x^2
(* b^2 A[x,x]+2 b c A[x, x^2]+c^2 A[x^2,x^2] *)

• Did you forgot the mixed term 2 b c A[x, x^2], in your example, or should it vanish for some reason? – jkuczm Oct 10 '15 at 13:05
• No, I forgot it. – simon Oct 10 '15 at 18:32
• When you posit that the operator is multilinear, that brings with it a notion of a base field (or at least a ring), that is, "scalars", and module elements. In the example provided, a,b,c are all scalars and x is a module element (e.g. a vector, or member of whatever space the operator works over). I bring this up because any implementation will need a way to distinguish one from the other. – Daniel Lichtblau Oct 11 '15 at 16:34

# Orderless

Orderlessness is easy just set Orderless attribute on your operator (I'll call it opA instead of A since symbol names starting with capital letters can collide with built-ins):

ClearAll[opA]
SetAttributes[opA, Orderless]


Now arguments of opA will be sorted to canonical order:

opA[b, a]
(* opA[a, b] *)


Linearity can be split into two parts: additivity and homogeneity.

Additivity is quite easy we need to add definition mapping our operator over expressions with head Plus:

opA[x_Plus, y_] := opA[#, y] & /@ x


Due to Orderless attribute, above definition of opA is tried with all permutations of arguments, so although we defined additivity only with respect to first argument, our operator is additive with respect to all arguments:

opA[a + b, c + d + e]
(* opA[a, c] + opA[a, d] + opA[a, e] + opA[b, c] + opA[b, d] + opA[b, e] *)


# Homogeneity

Homogeneity is a bit more difficult since we need to somehow tell Mathematica which expressions should be considered "scalars". This part depends strongly on your preferences and specific use case.

Let's define general function that answers question whether something is scalar, for purpose of linearity of our operator, or not.

ClearAll[scalarQ]
scalarQ[expr:_[__]] := FreeQ[expr, _?(Not@*scalarQ), {1}]
scalarQ[_] := True


scalarQ function when given expression with some arguments tests whether all arguments are scalars, when given any other expression it returns True.

By default all expressions are treated as scalars:

scalarQ /@ {1, a, a + b}
(* {True, True, True} *)


But we can tell Mathematica that certain symbol represents a non-scalar:

scalarQ[x] ^= False;


Now expressions containing this symbol will be considered non-scalars:

scalarQ /@ {x, x + a, 2 + a^x}
(* {False, False, False} *)


Depending on particular use case we may decide that some combinations of non-scalars give a scalar. For example Dot-products. By default they are considered non-scalars:

ClearAll[x, y]
(scalarQ[#] ^= False) & /@ {x, y};
scalarQ[x.y]
(* False *)


But we can add definition for Dot:

scalarQ[HoldPattern@Dot[__]] = True;
scalarQ[x.y]
(* True *)


Now let's go back to homogeneity of our operator and add definition for argument being multiplication of scalars by something else:

opA[Longest[scalar__?scalarQ] rest__, y_] := scalar opA[Times[rest], y]


Again due to Orderless attribute definition for first argument is enough.

ClearAll[a, b, c, x, y, z]
(scalarQ[#] ^= False) & /@ {x, y, z};
opA[a b x, c y z]
(* a b c opA[x, y z] *)


Notice that we declared x, y and z as non-scalars; a, b and c where treated as scalars.

# Putting it together

ClearAll[opA, scalarQ]

scalarQ[HoldPattern@Dot[__]] = True;
scalarQ[expr : _[__]] := FreeQ[expr, _?(Not@*scalarQ), {1}]
scalarQ[_] := True

SetAttributes[opA, Orderless]
opA[x_Plus, y_] := opA[#, y] & /@ x
opA[Longest[scalar__?scalarQ] rest__, y_] := scalar opA[Times[rest], y]


Example of usage:

ClearAll[a, b, c, x, y, z]
(scalarQ[#] ^= False) & /@ {x, y, z};
opA[x + 2 y + a b^5 z^2, x^2 + c y^3 z]
(* opA[x, x^2] + c opA[x, y^3 z] + 2 opA[x^2, y] + a b^5 opA[x^2, z^2] + 2 c opA[y, y^3 z] + a b^5 c opA[y^3 z, z^2]*)


and example from question:

opA[x, x] /. x -> 0 + b x + c x^2
(* b^2 opA[x, x] + 2 b c opA[x, x^2] + c^2 opA[x^2, x^2] *)