Fishing for monomials in a nested or partially factored polynomial stream

Question

I have a problem where I'd like to be able to take a multivariate polynomial whose variables are nonscalar and is not written explicitly as a sum of it's nonzero monomial terms and determine which monomials are present. For example,

z.(x+y.x)+y.(x.(y+1).z-x.z.y)

Since scalar multiplication is not present Coefficient, Expand and MonomialList are ineffective so a custom function is written based on Distribute to handle the expansion recursively through multiple depths. This is fine I suppose for the cases where I need all the monomials and their coefficients, but there are many times when I want to test for the presence of a few specific monomials, or to find the coefficients of those monomials.

Typically, the polynomials I produce of this type are of such a large degree or nested so deeply that it can take hours and gigs of memory just to be expanded out. I want to know if there is a better method for determining the coefficients of a given monomial without expanding the whole polynomial out. I haven't come across any literature that addresses such a problem, even in a more general computational setting. My best intuition so far is to traverse the expression like a tree, descending on branches which fit the criterion of the monomial and collecting coefficients along the way; unfortunately, implementing this in an efficient and "Mathematica-y" way has proven cumbersome.

---

If you want to run code through its paces, this will generate a polynomial not too dissimilar from my current applications

Nest[# + #.#.x + #.y.# &, z, 4]

You can try upping the recursion if you feel extra brave. If for whatever reason there's a desire to expand, the documentation for NonCommutativeMultiply includes an example of an expansion command which can be easily adapted to Dot. I do recommend adding the following rule otherwise it doesn't behave on sums of products.

ExpandNCM[h_Plus] := Map[ExpandNCM, h]

Edit: As requested, for wanting the coefficient of z.y.z.y.z.y.z.y.z from Nest[# + #.#.x + #.y.# &, z, 3], the following would give the result I desire

Clear[ExpandNCM]
ExpandNCM[h_Plus] := Map[ExpandNCM, h]
ExpandNCM[(h : Dot)[a___, b_Plus, c___]] := Distribute[h[a, b, c], Plus, h, Plus, ExpandNCM[h[##]] &]
ExpandNCM[(h : Dot)[a___, b_Times, c___]] := Most[b] ExpandNCM[h[a, Last[b], c]]
ExpandNCM[a_] := ExpandAll[a]

poly=Nest[# + #.#.x + #.y.# &, z, 3];
Coefficient[poly//ExpandNCM,z.y.z.y.z.y.z.y.z]

which is 10. This example takes relatively little time to compute, but if the Nest is bumped up to 4 iterations, the expansion takes substantially longer. The goal is to obtain the same result as Coefficient above without having to expand the full polynomial.

Answer 1

This might not be entirely correct or complete but it gives an idea of a recursive descent approach.

Clear[coeff];
coeff[dp_Plus, mon_, vars_] := coeff[#, mon, vars] & /@ dp
coeff[dp_, Dot[], vars_] /; FreeQ[dp, Dot] && FreeQ[dp, vars] := dp
coeff[a_.*dp_, mon_, vars_] /; MatchQ[dp, mon] && FreeQ[a, vars] := a
coeff[dp_, _, vars_] /; FreeQ[dp, Dot] := 0
coeff[a_*dp_Dot, mon_, vars_] /; FreeQ[a, vars] := 
 a*coeff[dp, mon, vars]
coeff[Dot[dp1_, dp2__], mon_Dot, vars_] := Module[
  {},
  Sum[coeff[dp1, Take[mon, j], vars]*
    coeff[Dot[dp2], Drop[mon, j], vars], {j, 0, Length[mon]}]
  ]
coeff[Dot[dp1_, dp2__], mon_, vars_] := 
 coeff[dp1, mon, vars]*coeff[Dot[dp2], 1, vars] + 
  coeff[dp1, 1, vars]*coeff[Dot[dp2], mon, vars]
coeff[dp_, x_, vars_] := Coefficient[dp, x]

Here is an example from the dot-based polynomial in the question.

vars = Alternatives @@ {y, z};
dotpoly = Nest[# + #.#.x + #.y.# &, z, 4];
monom = z.y.z.y.z;

coeff[dotpoly, monom, vars]

(* Out[317]= 12 *)