What is the most concise way to parallelize over all Table components at the same time?

Question

Question statement: Consider a Table with any number of indices {i,1,imax},{j,1,jmax},{k,1,kmax},.... Now I want to apply the same independent operation to each one of the elements while parallelizing the computation. What is the most concise way to do so?

---

Context:

So far, I have been using ParallelTable. However, this will separately take the lowest-order sub-array and paralellize over it separately. This may end up using fewer kernels at a time than available and slow down the process.

Consider the following example of a Table with 3 rows and 3 columns that I want to simplify in parallel:

ExTab = Table[..., {i,1,3},{j,1,3}];
SimplExTab = ParallelTable[FullSimplify[ExTab[[i,j]]],{i,1,3},{j,1,3}];

Now I have 10 kernels available, but ParallelTable will only use 3 kernels at a time, since it will only parallelize over each row at a time, even though it could use 9 and cut the time to roughly a third.

Of course, I could use Flatten on ExTab, apply ParallelTable, and reconstruct the structure of SimplExTab in the next step. This seems like an inelegant exercise, however. There should be a more concise way to do this. What is it?

---

---

More concrete example:

Let me give a little bit more precise example of my use case. Consider an example tensor

ExTens[n_,k_]:= Module[ {IndexStringList,TableString},
 IndexStringList = Table[ToString[{ToExpression["i:<>ToString[l]],1,k}],{l,1,n}];
 TableString = "Table[ExpandAll[Product[(x - RandomInteger[100])^RandomInteger[10],{m,1,18}]]";
 Do[TableString = TableString<>","<>IndexStringList[[l]],{l,1,n}];
 TableString = TableString<>"]";
 ToExpression[TableString]
]

For a given $n,k$, this function will randomly generate an $n$-index table (or tensor) in dimension $k$. Each component is a polynomial with integer coefficients, and variable powers of $x$ appearing. The thing that is not entirely obvious from the generated expressions is that every component can be factored using at most 18 positive integer roots of various multiplicity.

Now I want to recover the fact that the polynomials are obviously quite simple. Mathematica is able to do this by threading over the Table sequentially:

RandomTens = ExTens[5,5]
FullSimplify[RandomTens]

This does give the simple result of each component factorized by using the 18 integer roots. However, you start to notice the non-zero computation time even in this reasonably simple case ($\sim$ minutes on my machine).

---

Attempt at fully parallelizing concrete example:

I attempted to accelerate the example by using ParallelMap[] as suggested by Daniel Huber in the comments. The idea is that one can use ParallelMap to apply a function to each component of a simple Table:

ParallelMap[f,{{a,b},{c,d}},{-1}];
-> {{f[a],f[b]},{f[c],f[d]}

This takes all parts of the table {{a,b},{c,d}} of depth one (which in this case are the table components), and applies $f$ to them. However, this stops working for RandomTens, since its components are expressions that have a longer depth due to being polynomials. So I instead use:

 ParallelMap[FullSimplify[#]&,RandomTens,{ArrayDepth[RandomTens]}];

This does parallelize the computation partially, but surprisingly enough, it is only as good as the ParallelTable[] method given above! In particular, it never uses the full 10 cores, it launches only $k$ cores, and goes again "row by row" in some sense!

Answer 1

After exploring a few options in Mathematica and elsewhere, I came to the conclusion that there is probably no in-built function or concise expression that would do what I am looking for (speaking of v13.2 MMa here).

I believe that the reason for this is that a general array or list can have a rather complex and varied structure, and brute-forcing the parallelization over all components may or may not give a better result than the default, conservative Mathematica treatment.

I have written my own function that does what I mention in the OP:

ParTabSimp[f_,T_] :=   Block[{FlatT,Tind,Tlength,SubT,Tcomps,FlatSimpT,k,IndexList,IndexSum},
  FlatT = Flatten[T];   (*Flattened version of T*)
  Tind = ArrayDepth[T]; (*Number of indices of tensor*)
  SubT = T;             
  Do[Tlength[i] = Length[SubT]; SubT = SubT[[1]],{i,1,Tind}]; 
  (*The assumption is that T has a tensor-like structure, that is, it is a table of i1*i2*... numbers, no variable-length structures etc. are assumed*)
  Tcomps = Product[Tlength[i],{i,1,Tind}]; (*Tlength[i] tells you the range of the i-th index*)
  FlatSimpT = ParallelTable[f[FlatT[[i]]],{i,1,Tcomps}];
  IndexSum = " ";IndexList = " ";
  Do[
    IndexSum = IndexSum<>"+i"<>ToString[j]; 
    IndexList = IndexList<>",{i"<>ToString[j]<>",1,"<>ToString[Tlength[j]]<>"}"
  ,{j,1,Tind}];
  ToExpression[("Table[FlatSimpT[["<>IndexSum<>"]]"<>IndexList<>"]")]
]

This function breaks up a multi-dimensional table T into its single entries, and applies the function f[] on each component while using all kernels available, then reconstructs the original tensor structure.

The use of the function is for example:ParTabMap[FullSimplify[#]&, RandomTens], using the random tensor RandomTens from the OP, and this can provide a significant edge over ParallelMap[FullSimplify[#]&,RandomTens,{ArrayDepth[RandomTens]}] in some use cases.