Efficient way to apply linear function to multivariate polynomial

Question

Suppose I start with an expression that is a multivariate polynomial in $x_k$'s, $$W = a + b \cdot x_1^{n_1} x_2^{n_2} x_3^{n_3} x_4^{n_4} + c \cdot x_1^{m_1} x_2^{m_2} x_3^{m_3} x_4^{m_4}$$

where $a, b, c$ are just constants, and $n_i, m_j$ are nonnegative integers, and $x_k$'s are real scalars.

Now, suppose I apply a linear operator $L$ (what I have in mind is the expectation operator from probability and statistics) to the above, and associate $$L[ x_1^{n_1} x_2^{n_2} x_3^{n_3} x_4^{n_4} ] = M( n_1, n_2, n_3, n_4)$$

(i.e. in probability and statistics, $M$ would correspond to the moment $E[X_1^{n_1} X_2^{n_2} X_3^{n_3} X_4^{n_4}]$), so that $$L[W] = a + b \cdot L[x_1^{n_1} x_2^{n_2} x_3^{n_3} x_4^{n_4}] + c \cdot L[x_1^{m_1} x_2^{m_2} x_3^{m_3} x_4^{m_4}] \\= a + b M(n_1, n_2, n_3, n_4) + c M(m_1, m_2, m_3, m_4)$$,

and suppose that the function $M$ can be computed numerically quickly. The objective is to compute $L[W]$ efficiently.

Question: Wow does one do this efficiently when there are many, many terms in the polynomial?

That is, suppose

w = a + b * x1^n1 * x2^n2 * x3^n3 * x4^n4 + c * x1^m1 * x2^m2 * x3^m3 * x4^m4 

The current method that I have is effectively to use CoefficientRules applied to $W$, and say get wcr = CoefficientRules[w], and then compute Keys[wcr], which will then get me a list of the exponents associated with $x_k$'s. Apply my function $M$, say mlist = mfun @@@ Keys[wcr]. And then finally, to put everything back together, I will then simply multiply, and get resultL = Values[wcr] * mlist (I could also use Dot, I suppose). But in the application I have in my mind, I have potentially hundreds and thousands of terms, so I figured out that the multiplication step is taking considerable amount of time.

Answer 1

One obvious answer is to use the answer to your related earlier question with the straightforward modification of renaming the random variables:

Clear[expect]

expect[expr_Plus] := Map[expect, expr]

expect[Times[x_, y__]] /; (FreeQ[x, x1 | x2 | x3 | x4]) := 
 x expect[Times[y]]

expect[Times[x_expect, y__]] := x expect[Times[y]]

expect[expr_?(FreeQ[#, x1 | x2 | x3 | x4] &)] := expr

expect[Power[expr_Plus, n_]] := expect[Expand[Power[expr, n]]]

expect[Power[x_expect, n_]] := x^n

The apply it as follows:

w = a + b*x1^n1*x2^n2*x3^n3*x4^n4 + c*x1^m1*x2^m2*x3^m3*x4^m4

expect[w]

(*
==> a + c expect[x1^m1 x2^m2 x3^m3 x4^m4] + 
 b expect[x1^n1 x2^n2 x3^n3 x4^n4]
*)

Now you just have to add a definition for the expectation value of the monomials as follows:

expect[x_] := M[x]

expect[w]

(*
==> a + c M[x1^m1 x2^m2 x3^m3 x4^m4] + 
 b M[x1^n1 x2^n2 x3^n3 x4^n4]
*)

Here, I assume the fast implementation of the expectation value you mention in the question is called M. By omitting any restrictions on the pattern in expect[w_], this definition is the last one that kicks in, after the linearity properties have already been applied.