Parallel Experiments

Roman E. Maeder

MathConsult Dr. R. Mäder
Samstagernstraße 58a
8832 Wollerau
Switzerland
http://www.mathconsult.ch

Abstract

The developer of Parallel Computing Toolkit presents examples of large-scale parallel calculations and discusses performance and reliability issues. Ongoing experiments with irreducibility proofs of polynomials with a total CPU time of several years were used to develop parallel search methods and tools for unattended computations.

© 2005 MathConsult Dr. R. Mäder. All rights reserved.

Agenda

Infrastructure, with Special Thanks to Orion Multisystems

({{, 
In this demonstration, the frontend runs on my laptop computer,
the master kernel and slave kernels run on an Orion Multisystems DT-12,
a desktop cluster with 12 processors.
Orion made such a machine available to me shortly after last year's conference.}})

The results presented here are the result of almost one year of computation, or equivalently 12 CPU years, with PCT V2 and Mathematica 5.0—5.2.

\(TableForm[
      RemoteEvaluate[{$ProcessorID, $MachineName, $SystemID, $Version, 
          TimeUsed[], \ MaxMemoryUsed[]}], 
      TableHeadings -> {None, {"ID", "host", \ "OS", "Mathematica Version", \ "CPU Time", "Memory"}}]\)

The Problem: Irreducibility of a Family of Polynomials

Are the polynomials (∑_(i = 1)^(n + 1))i (x^(i - 1)) all irreducible? Here we use the equivalent reciprocal form (n + 1) + 
        n x + (n - 1) (x^2) + … + 2  (x^(n - 1)) + (x^n).

poly[n_, x_] := (∑_(i = 1)^n)(n - i + 1) (x^i) + n + 
        1

p10 = poly[10, x]

In the following we explore various ways to search for counterexamples, that is, nontrivial factorizations.

This question was posed by Michael Filaseta in 1986.

Interactive Search for Factors

We can use FactorList, which gives a list of factors and their multiplicity.

FactorList[p10]

The number of factors is equal to the sum of the exponents, disregarding the numerical content (which is always 1 here).

nFactors[poly_]:=Plus@@Rest[FactorList[poly][[All,2]]]

nFactors[p10]

A Simple Dispatcher

You make the definitions known to all remote kernels.

ExportEnvironment[{poly, nFactors}];

Now you can start the computation, requiring that it print each result as it continues. To stop the computation, abort it by choosing Kernel ⊳ Abort Evaluation or pressing ⌘[.].

d = 500;  ids = {}; results = {}; 
While[True,
While[$QueueLength <= 0, 
AppendTo[
                ids, Queue[{d}, {d, 
                    nFactors[
                      poly[d, 
                        x]]}]]; 
d++; 
If[! QueueRun[],  Break[]]; 
];
If[
              Length[ids] > 0, 
{res, id, ids} = 
                WaitOne[ids]; 
Print[
                res]; 
AppendTo[results, 
                res]; 
];
]

.

We need to clean up after an interrupt.

ResetSlaves[];

Handling the Abort

Make sure the ResetSlaves[] is not forgotten.

d = 500;  ids = {}; results = {}; 
CheckAbort[
\(While[
            True, 
While[$QueueLength <= 
                0, 
AppendTo[ids, 
                Queue[{d}, {d, 
                    nFactors[
                      poly[d, 
                        x]]}]]; 
\(d++\); 
If[\(! QueueRun[]\), \ 
                Break[]]; 
]; 
If[
              Length[ids] > 0, 
{res, id, ids} = 
                WaitOne[ids]; 
Print[
                res]; 
AppendTo[results, 
                res]; 
]; 
]\),
  
\(Print["resetting..."]; 
ResetSlaves[]\)
];
Length[
      res]

Refinements

We can easily exclude primes that give linear factor by trial division by all elements of GF[p], or equivalently by evalution.

With[{p = 11, 
        n = 10}, 
Mod[poly[n, #] & /@ Range[0, p - 1], 
        p]
]

Excluding higher-degree factors is a bit more complicated and will be discussed later.

Batch Processing

This approach is already described in an example that is part of the Parallel Computing Toolkit documentation.

Save Results

This is easily done with PutAppend[].

…
{6760.960178, 3017, 1, 5}
{8334.958893, 3014, 1, 2}
{8387.436915, 3016, 1, 4}
{9438.018203, 3015, 1, 3}
{10487.562649, 3019, 1, 7}
{10889.343568, 3021, 1, 9}
…

Print regular progress reports

A number of Print[] statements and a global variable that determines the amount of information to print.

Degree 6792 queued, 12 results received in 1426m, 63247s avg (mex: 6762).
Degree 6804 queued, 12 results received in 1040m, 82574s avg (mex: 6776).
Degree 6816 queued, 12 results received in 773m, 54238s avg (mex: 6776).
Degree 6828 queued, 12 results received in 1097m, 50822s avg (mex: 6804).
…

Restarting

To restart the computation after a power failure, read the previous results and figure out where to continue.

Results (up to degree 8253)

The full power of Mathematica is available to analyze the results received.

res = 
      ReadList["/home/orion/clusterfs/maeder/results.m"]; 
Length[res]

Graph of runtimes vs. degree.

Obviously this graph suggests certain possible tweaks in the factorization algorithm.

Hard cases by runtime

\(TableForm[\({Round[#1/3600*10]/10.0, ##2} &\) @@@ 
        Take[Reverse[Sort[res]], 12], 
      TableHeadings -> {Automatic, {"time (h)", "degree", "factors", "$PID"}}]\)

Easy cases (among last 10%)

\(TableForm[Take[Sort[Take[res, \(-Floor[Length[res]/10]\)]], 10], 
      TableHeadings -> {None, {"time (s)", "degree", "factors", "$PID"}}]\)

Total computation time

Plus @@ degtime[[All, 2]]

Days

%/(24*60*60)

Approx. real time in days

%/12

Alternative: Modular Factorizations

A full factorization does too much work. It switches to an expensive integer factor reconstruction after trying a small number of primes. In most cases we can prove irreducibility from a number of modular factorizations alone.

Factor Degree Sets

The degrees of all irreducible factors mod Prime[pi].

factorDegrees[d_, pi_]:=
Module[{t, faclist, x, degs},
faclist=FactorList[poly[d, x],Modulus -> Prime[pi]];
degs=Flatten[Table[Exponent[#1, x],{#2}]&@@@faclist];
degs
]

Let's try it

d = 1180;

d1 = factorDegrees[d, 1]

All possible degrees

Set of all possible degrees, all possible subset sums (but computed iteratively).

degreeSet[degs_]:=Fold[Union[#1,#1 + #2]&,{0},degs]

s = s1 = degreeSet[d1]

Next attempt

d2 = factorDegrees[d, 2]

s2 = degreeSet[d2]

s = Intersection[s, s2]

Next one

d3 = factorDegrees[d, 3]

s3 = degreeSet[d3]

s = Intersection[s, s3]

QED

Avoiding bad primes

Omit primes with factors of low degrees in the hope of obtaining fewer factors and smaller degree sets.

Testing for factors of degree k: evaluate the polynomial at all points of GF[(p^k)]. If at least one result is 0, the polynomial has a factor of degree k (or one of k's factors).

I'm Feeling Lucky

Luckily the Horner formula gives all intermediate polynomials

(… ((0  x + 1) x + 2) x + … + n) 
        x + (n + 1)

Rest[FoldList[Expand[#1*x + #2] &,0,Range[5]]]

Therefore, we can compute the evaluations incrementally.

Evaluation

Evaluate poly[d,x] mod p at a

With[{d = 10, p = Prime[5], a = 1},
Rest[FoldList[Mod[#1*a+#2,p]&,0,Range[d + 1]]]
]

All points in one go

a can be a list of points, the code is the same.

With[{d = 10, p = Prime[5]},
With[{a = Range[0, p - 1]},
Rest[
              FoldList[Mod[#1*a + #2, p] &, 0, Range[d + 1]]]
]]

Which lists contain a zero?

Times @@@ %

All degrees that have a linear factor:

Flatten[Position[%, 0]]

Higher-degree factors

Arithmetic in GF[(p^k)]: represent elements as lists of length k, considered coefficients of a primitive polynomial.

Compute the coefficients of x*a + i symbolically beforehand and compile the result. Note that a is the same everytime, and be inserted right into the expression.

Selection of extension degrees

As the number of elements of GF[(p^k)] grows quickly with k, rather than use the lcm of the desired degrees (what I am doing now), we could compute for several values of k.

For example, to cover all degrees up to 10, we need these extensions: 10,9,8,7,6. All lower ones are divisors of these.

Another consideration is that factorization times grow with p, so we should not be too strict. For p>17, I check only for linear factors.

Data Types

Using data types in the way I advocate in my books helps keep the code well organized.

Modular Factorizations

Keep track of modular factorization results. Attributes include the degree, prime, factor degree set, and time taken.

new[modular][d_, pi_, degrees_, 
        time_] := 
modular[d, pi, degrees, time]

Selectors

modular /: deg[modular[d_, pi_, degrees_, time_]] := d
modular /: pi[modular[d_, pi_, degrees_, time_]] := pi
modular /: degrees[modular[d_, prime_, degrees_, time_]] := 
      degrees
modular /: time[modular[d_, prime_, degrees_, time_]] := 
      time

Auxiliary functions

modular /: degreeSet[m_modular] := 
      Fold[Union[#1, #1 + #2] &, {0}, degrees[m]]

Polynomial Factorizations

Keep track of the progress for one degree. Attributes include the list of modular factorizations and the surviving degree set.

new[factorization][
        d_Integer?Positive] := 
Module[{modulars = {}, 
          degset = Range[0, d]}, 
factorization[d, 
          modulars, degset]
]

Selectors

factorization /: degree[factorization[d_, modulars_, degset_]] := 
      d
factorization /: modulars[factorization[d_, modulars_, degset_]] := 
      modulars
factorization /: degset[factorization[d_, modulars_, degset_]] := 
      degset

Auxiliary functions

factorization/:addModular[factorization[d_, modulars_, degset_], m_modular]:=(AppendTo[modulars, m];
degset=Intersection[degset, degreeSet[m]];
)

factorization /: irreducibleQ[f_factorization] := 
      Length[degset[f]] ≤ 2

Database

We simply store a factorization object for each degree as a downvalue.

We can use dynamic programming for the initialization of the store as needed (when one degree is access for the first time).

fstore[d_Integer]:=f  store[d] = new[factorization][d]

Whenever a new result arrives from one of the remote kernels, the corresponding fstore[d] element is looked up and augmented with the new information.

The Dispatcher

The main loop running on the master kernel.

Progress Report

3321/7 received, 183s, 14 factors, sofar 1 primes, 3252 degrees remaining
GF[5^6] extended to 3423, 13s, 1 polynomials released
GF[23^1] extended to 3421, 0s, 1 polynomials released
GF[11^4] extended to 3423, 7s, 1 polynomials released
3315/10 received, 499s, 7 factors, sofar 3 primes, 4 degrees remaining
3320/3 received, 247s, 6 factors, sofar 2 primes, 12 degrees remaining
3313/13 received, 379s, 8 factors, sofar 4 primes, 0 degrees remaining
3313 IRREDUCIBLE, 4 primes, {8, 10, 11, 13}, 1147s total
…

Results

Runtime vs. degree (with a quadratic fit)

Breakdown by number of moduli

Interesting Cases

Data mining in the area of modular factorizations.

Hard Cases, by runtime

Easiest among last 100

Lucky ones, only one trial

(There are more of those in the degree range from 2000—3000)

Stubborn Ones

Modular Factorizations

Lucky ones, irreducible

Unlucky ones

many factors,excluding Prime[1]==2

Hot off the CPUs

Batch Tools

Collect all code in a notebook. The initialization cells contain

Other sections may contain

The Main Loop

The main loop is not run automatically. It is only set up.

mainloop := 
Module[
  
]

In this way evaluating the initialization does not start the computation (important if the notebook is used to analyze results only).

The Package

Each time the notebook is saved, all initialization cells are written to a package .m file.

The batch script

A Perl script starts a Mathematica kernel and feeds it the necessary commands.

Mathematica output is redirected to a file. This is where progress reports will be collected.

Starting a batch computation

The computation is started from a shell. The Mathematica command inserted reads the package and then invokes the main loop.

echo 'Get["PolyModular.m"];mainloop' | batchrun -batch - >&batch.err &

(Error) output from the script is redirected to a file.

Stopping a batch computation

In each trip through the main loop the program checks for the presence of a particular file. If it is there, the main loop exits.

Reliability

The remote kernels will grow in size over time. If one exceeds the virtual memory limit, it will quit and the failure recovery of PCT will assign its task to a different kernel. After this happens, I usually stopped the computation and started over.

Conclusions

Resources on the Web

Having a personal cluster on one's desktop can be a lot of fun.

The experiments gave me plenty of ideas for new features of the next version of the Parallel Computing Toolkit.

CloseSlaves[]

Quit[]