Why is Flattening a CoefficientArray so slow?

I have a vector of $n$ degree $n$ polynomials in $(x,y,z)$, each of whose coefficients is an inhomogenous linear expression in $3 n^2$ variables. I want to write down the linear equations in these $3 n^2$ variables which say that all these coefficients are identically zero. Here $n$ is about $15$. I used to be doing this with a combination of various Map[], Coefficient[] and CoefficientList[] commands, but the result was much slower than I wanted.

Tonight, I discovered the CoefficientArrays[] command in the documentation, which seemed built for exactly this job. To my surprise, I found that replacing one of my Coefficient[] calls with CoefficientArray[] was a dramatic improvement, but replacing the other CoefficientList[] was a disaster.

Of course, I know what to do: Replace the good one but not the bad one! But I am curious why.

The following code sets up test data. In the actual use, f and PolyBasis are not random polynomials but part of a larger computation, but I found that I could use random data for profiling without changing the speed:

n = 15;

f = Sum[Random[] x^i y^j z^(n - i - j), {i, 0, n}, {j, 0, n - i}];
PolyBasis = Table[Sum[Random[] x^i y^j z^(n-1 - i - j), {i, 0, n-1}, {j, 0,
n-1 - i}], {n}];

AA = Table[a[i, j], {i, 1, n}, {j, 1, n}];
BB = Table[b[i, j], {i, 1, n}, {j, 1, n}];
CC = Table[c[i, j], {i, 1, n}, {j, 1, n}];

MM = x AA + y BB + z CC;

LeftKer = MM.PolyBasis - Prepend[ConstantArray[0, n - 1], f];


Here is what I used to be doing:

AbsoluteTiming[EqL = Flatten[Map[CoefficientList[#, {x, y, z}] &, LeftKer]];]
(* {0.562461, Null} Kind of slow. *)

shortEqL = Select[EqL, ! (# === 0) &];

vars = Flatten[Join[
Table[a[i, j], {i, 1, n}, {j, 1, n}],
Table[b[i, j], {i, 1, n}, {j, 1, n}],
Table[c[i, j], {i, 1, n}, {j, 1, n}]]];

AbsoluteTiming[Big = Map[Coefficient[#, vars] &, shortEqL];]
(* {20.197230, Null} REALLY SLOW *)

AbsoluteTiming[small = shortEqL /. Map[(# -> 0) &, vars];]
(* {3.257321, Null} Significantly slow *)


Replacing the last two commands with CoefficientArrays[] is a huge gain:

AbsoluteTiming[{small2, Big2} = CoefficientArrays[shortEqL, vars];]
(* {0.105790, Null} Yippee! *)


Replacing the "kind of slow" command with CoefficientArrays gains a little:

AbsoluteTiming[EqL2=CoefficientArrays[LeftKer, {x, y, z}];]
(* {0.401718, Null} *)


But Flatten[EqL2] runs for many minutes without stopping, as does Normal[EqL2]!

Question: Why is Flatten[]ing the output of CoefficientArrays[] so bad? And is there something smart I should be doing to improve on EqL = Flatten[Map[CoefficientList[#, {x, y, z}] &, LeftKer]];?

For the curious, we are trying to remove some of the bottlenecks in the algorithm described here.

-
The link to math.uni-konstanz.de/~plaumann/PSSV12.pdf seems to be broken. –  A.G. Dec 18 '13 at 15:46

When you are flattening the expression, you are asking Mathematica to reserve memory for at least

Total@Cases[EqL2, SparseArray[x_, y_, z__] :> Times @@ y]
(*322850400 elements *)


To be more precise, you are only filling

Length@Select[Flatten@Cases[EqL2, s_SparseArray :> ArrayRules@s] /.
Rule -> List, ! FreeQ[#, a | b | c] &]
(* 2040 elements*)


with your coefficients, and the rest seems being filled with zeroes.

Perhaps you may use ArrayRules[] and keep the Sparse Arrays instead of normalizing them.

-