Extracting coefficients of very large polynomials (perhaps by improving Expand or ExpandAll)

Question

Is there a way to improve or bypass Expand (or ExpandAll) for extremely large polynomials?

I have a polynomial system of the form $$ 0 = \sum_n a_n d_{i_1,i_2} \cdots d_{i_{p-1},i_p},$$ in which

• the $d_{i_k,i_j}$'s are fixed
• $i_k \in \{ 1,2,3,4,5,6 \}$
• the $a_n$ are undetermined.

I solve for the $a_n$ by collecting the coefficients of the $d$'s, and solving the resulting system of equations among the $a$'s.

My problem occurs with the very large system modeled as:

(* Argument of the d's *)
dr := RandomInteger[{1, 4}];
(* Label the a's *)
ar := RandomInteger[{1, 20}];
(* Numerical factors may show up *)
cr := RandomInteger[{-10000, 10000}];
pr := RandomInteger[{-4, 4}];
expr =
Sum[
    Sum[cr*a[ar]*d[dr, dr] d[dr, dr], {20}]*
    Product[Sum[pr*d[dr, dr] d[dr, dr], {4}], {13}]
    ,
    {15}
    ];

My system always comes from combining a sum of rational terms, and taking just the numerator, so they always have the form (modeled by the code):

$$ 0 = \sum_{i} \textrm{numerator}_i \times \prod_{j \ne i}\textrm{denominator}_j, $$

where the numerators and denominators are all polynomials in the $d$'s, and the $a$'s appear in the numerators. I mention this because perhaps it can be used in bypassing my problem.

My problem is in

expr = ExpandAll[expr];

or

expr = Expand[expr];

for expressions that are this long. This is a huge time and memory drain (which is totally fair, but it's also the problem I'm hoping to bypass).

I have no problem in constructing the coefficient equations after this (using GatherBy), and one of the reasons I need to expand is because the GatherBy equation construction is incredibly efficient.

I'm happy to use a different method to construct the equations if it will drastically reduce the net time and/or memory.