Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

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

1 Answer 1

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.