# Fishing for monomials in a nested or partially factored polynomial stream

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.

• Could you please post an usage example using Nest[# + #.#.x + #.y.# &, z, 4] and the desired result? Jul 28, 2015 at 22:09
• You meant Times and not Dot, right? I ask because otherwise it's not a polynomial in the usual sense. Jul 28, 2015 at 22:43
• @DanielLichtblau Nope. The provided examples use Dot. Jul 28, 2015 at 22:46
• Sooo...what are the "monomials"? Things like x.y.x.x.z.y with order being important? Jul 28, 2015 at 22:52
• @DanielLichtblau Correct. The variables are noncommutative in multiplication. Jul 28, 2015 at 23:02

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 *)

• This is an interesting approach. I can't find any case so far where it fails and it's much faster than expanding for large problem sizes (even comparable in small cases). I'll test it more rigorously but for now you have my appreciation and my +1. Aug 5, 2015 at 4:02